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This is a final report that summarizes the results achieved under this grant The first 
j major accomplishment is the development of the sublaminate modeling approach and shear 
r deformation theory. The sublaminate approach allows the flexibility of considering one ply 
jj or groups of plies as a single laminated unit with effective properties. This approach is 
} valid when the characteristic length of the response is small compared to the sublaminate 
| thickness. The sublaminate approach was validated comparing its predictions with a finite 
i element solution [J], A shear deformation theory represents an optimum compromise 
! between accuracy and computational effort in delamination analysis of laminated 
I composites [2]. This conclusion was reached by applying several theories with increasing 
| level of complexity to the prediction of interlaminar stresses and strain energy release rate in 
1 a double cracked-lap-shear configuration. 

The shear deformation theory and sublaminate approach was applied to the free- 
edge delamination[I,i] and internal delamination analysis [4\ of laminated plates including 
the influence of hygrothermal stresses [5,6] and combined loading [7]. the analysis was 
also applied to tapered laminates subjected to tensile loading [8,9]. 

The second accomplishment is the development of the variationally asymptotical 
analysis for thin-walled anisotropic beams with closed cross sections [ 10-12 ] . The theory 
is a prerequisite for isolating the influence of damage by comparing predictions with an 
reference undamaged configuration. Existing composite beam theories have significant 
differences in the derived expressions for the stiffness coefficients. The variationally 
asymptotical analysis was developed in order to isolate the effects contributing to these 
differences. The major advantage of this approach lies in the fact that the displacement field 
is not assumed a priori as is the case for the existing theories and emerges as a result of the 
analysis. Moreover, the assumed displacement fields in the existing theories follow the 
classical isotropic formulation. However, no proof is provided with regard to the validity 
of such a displacement field for anisotropic materials. 

The displacement field which resulted from the theory showed two new 
contributions which were identified as out-of-plane warping due to axial strain and 
bending. These contributions emerge in addition to the classical out-of-plane torsional 
warping and are significantly influenced by the material's anisotropy. However, they 
vanish for materials that are orthotropic or whose properties are antisymmetric relative to 
the beam middle surface. These configurations coincide with the cases where the 
predictions of the existing theories are in agreement with test results and numerical 
simulations. For generally anisotropic materials the error associated with the existing 
theory predictions correspond to the neglect axial strain and bending related out-of-plane 
warping. 

In addition to providing a definitive answer to the reasons for the disparity in 
existing theories predictions, the variationally aymptotical theory provides a consistent 
approach to deriving the displacement field in anisotropic structures. A number of 
investigators have now adopted this approach for the modeling of initially curved and 
twisted composite beams and laminated composite plates[13, 14]. Moreover, the closed 


form expressions indicate that the new contributions are proportional to the extensional 
strain and bending curvature. This provides a proof for the work of Kosmatka [15] where 
an improvement to the displacement field was proposed by adding two terms which are 
proportional to the extensional strain and bending curvatures. However, their contributions 
were determined using a finite element simulation. 

The details of the sublaminate and Variationally asymptotical analyses are provided 
by the work of Ref. 12 which is provided in Appendix A for convenience. A list of the 
publications and presentations related to the Grant is provided in Appendix B. 
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CHAPTER I 
INTRODUCTION 


1.1 Background 

The use of fiber reinforced composites is increasing in engineering applications. One 
of the major issues in composite structures is the understanding of the role of the ma- 
terial’s anisotropy on the deformation modes, damage modes and failure mechanisms. 
This research work addresses these stiffness and strength related issues by developing 
analytical models for the prediction of deformation modes and their coupling effects 
and damage onset and growth in laminated composites. Accurate prediction of stiff- 
ness, response, damage modes and failure mechanisms is bound to lead to the design 
of efficient and damage tolerant composite structures. 

Delamination is a predominant failure mode in continuous fiber reinforced lam- 
inated composite structures. Based on the location and direction of growth, there 
are two distinct types of delamination, namely, free edge delamination and local or 
transverse crack tip delamination. In many cases, both types occur concurrently with 
varying levels of interaction. 

In the first part, of this work shear deformation models including hygrothermal 
effects are developed for the analysis of mid-plane edge delamination and local de- 
lamination originating from transverse cracks in 90° plies. The results of these models 
are combined with a previously developed shear deformation model for mixed-mode 
edge delamination to yield a unified analysis of delamination and the ability to iden- 
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tify the critical failure modes and loads. 

Elastically tailored composite design are being used to achieve favorable defor- 
mation modes under a given loading environment. Coupling between deformation 
modes such as extension-twist, or bending-twist is created by an appropriate selec- 
-tion of fiber orientation, stacking sequence and materials. An example is the X-29 
swept forward wing aircraft where a laminated composite skin is used to create the 
bending-twist coupling required to handle divergence. This design uses AS-1/3501-5A 
graphite/epoxy wing covers with —45° outboard plies 9° forward of the wing’s 40 /( 
chord line. Elastically tailored composite rotor blades have the potential to be used 
in rotorcraft structures in order to control flapping and twisting motions at different 
rotor speeds. This concept can be utilized in tilt rotor aircraft in order to achieve a , 
compromise between hover performance and forward flight propulsive efficiency. A~- 
change in the blade twist between flight modes can be developed through the use 
of extension-twist coupling. as implemented in the XV-15 tilt rotor aircraft. Twist 
control is achieved by assuming a 15 percent change in operating rpm between hover 
and forward flight regimes. 

The coupling of deformation modes provides a flexibility to meet design require- 
ments on the aeroelastic behavior, dynamic response and stability of structures and 
results in improved fatigue life and durability. 

A prerequisite for the implementation of an elastically tailored concept, is the 
development of an analytical model which accurately predicts the various stiffness 
components and isolate the material and geometrical parameters controlling the be- 
havior. 

In the second part, a variationally and asymptotically consistent theory for thin- 
walled beams that incorporates the anisotropy associated with laminated composites 
is developed. The theory is based on an asymptotical analysis of 2D shell energy. 
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The major advantage of this approach lies in the fact that the displacement function 
is not assumed a priori and is determined as a result of the minimization of the energy 
functional. As a result, two non classical contributions to the warping emerge. While 
these new contributions vanish for isotropic and orthotropic materials, they have a 
significant influence on the response of generally anisotropic materials. The accuracy 
of previously developed theories is assessed by comparing the resulting displacement 
fields and an assessment of the significance of shear deformation is presented. Com- 
parison of predictions with finite element simulation and test results illustrate the 
consistency and accuracy of the developed theory. 

The delamination analysis model is present ed in the first part of this work, this is 
followed by the development of the thin-walled anisotropic beam theory. Each part 
includes a literature survey in order to place the present work in proper prospective. 
A comparison of prediction is presented in order to validate the developed theories 
and assess their accuracy. 


4 


CHAPTER II 

DELAMINATION ANALYSIS 

This chapter addresses damage modeling in laminated composite plates. A review 
of previous work is presented first, this is followed by a development of the analytical 
model. 

2.1 Review of Previous Work 

Failure in laminated composite materials often initiates in the form of matrix frac- 
tures, namely, transverse matrix cracks and delaminations. Based on the location and 
direction of growth, two distinct types of delamination can be discerned. These two 
types are called edge delamination and local or transverse crack tip delamination, as 
shown in Fig. 2.1. Edge delaminations initiate at the load free edges of the laminate 
whereas local delaminations start from a transverse matrix crack. Transverse matrix 
cracks refer to intralaminar failures whereas delaminations refer to interlaminar fail- 
ures. Transverse cracks usually occur within laminates where the fibers run at an 
angle to the primary load direction and hence the name. In many cases, both types 
occur concurrently with varying levels of interaction. 

It has been observed [1] in simple tension tests of uniform rectangular cross section 
specimen (Edge Delamination tests) that delaminations initiate along the load free 
edges and propagate normal to the load direction as shown in Fig. 2.1. Transverse 
matrix cracks running parallel to the fibers have also been observed in off-axis and 
90° plies. Such transverse cracks extend through the thickness of similarly oriented 
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plies and terminate where the ply orientation changes. Delaminations can also origi- 
nate at the interfaces where transverse cracks terminate. These transverse crack tip 
delaminations or local delaminations, grow normal to the transverse crack from which 
they originate. In the case of 90° plies, the growth direction is parallel to the load. 

The growth process of edge delaminations and local delaminations is often mod- 
eled using a fracture mechanics approach leading to the calculation of a strain energy 
release rate. This is because the strain energy release rate can correlate delamination 
behavior from different loading conditions and can account for geometric dependen- 
cies. The strain energy release rate associated with a particular growth configuration 
is a measure of the driving force behind that failure mode. In combination with ap- 
propriate failure criteria, the strain energy release rate provides a means of predicting 
the failure loads of the structure. 

Several methods are available in the literature for analyzing edge delaminations. 
These include finite element modeling as in J2], J3], and [4], the complex variable stress 
potential approach [5], a simple technique based on classical laminate theory [1] and 
a higher order laminate theory including shear deformations [6]. Finite element mod- 
els provide accurate solutions but involve intensive computational effort. Classical 
laminate theory (CLT) provides simple closed form solutions and is thus well suited 
for preliminary design evaluation. However, CLT provides only the total energy re- 
lease rate, and thus, in a mixed mode situation, there is insufficient information to 
completely assess the delamination growth tendency. A higher order laminate theory 
including shear deformations has the ability to provide the individual contributions 
of the three fracture modes while retaining the simplicity of a closed form solution. 
A shear deformation model is available for off-mid-plane edge delamination and has 
been shown to agree well with finite element predictions [7]. 

Crossman and Wang [8] have tested T300/93 4 graphite/epoxy [±25/90 n ], sped- 
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mens in simple tension and reported a range of behavior including transverse cracking, 
edge delamination and local delamination. O’Brien |9] has presented classical lam- 
inate theory solutions for these specimens, demonstrating reasonable agreement in 
the case of edge delamination but with some discrepancies in the local delamina- 
tion predictions. The local delamination model overestimates the failure strains for 
[±25/90,,], specimens for small values of n mainly due to the implicit critical strain 
energy matching used. 

A finite element model combining edge and local delaniinations has been pro- 
posed by Law [10]. His predictions, however, do not fully explain the dependency of 
the critical strain on the number of 90° plies. A similar three-dimensional finite ele- 
ment analysis including hygrothermal effects has been performed by Wang et al. [11] 
to determine the delamination onset load for combined delamination, qualitatively 
demonstrating stable crack growth. 

A three-dimensional finite element analysis of delamination from matrix cracks 
has been developed by Fish and 0’Brien[l2]. They conducted an experimental and 
analytical study on the influence of matrix cracking on delamination in [±15/ — 

' 90„/ — 15], glass-epoxy laminates subjected to monotonically increasing tension loads. 
Experimental results showed that local delaminations form at the intersection of 
matrix cracks in the ±15° plies and the free-edge. Comparison of a Quasi-three- 
dimensional (Q3D) finite element results with a three-dimensional (3D) finite element 
analysis showed significant differences in the relative and absolute magnitudes of the 
interlaminar stress components. Thus, discrepancies in failure predictions may exist 
between Q3D and 3D analysis. The results of this study emphasized the importance 
of incorporating the various damage mechanisms that influence subsequent damage 
development in the failure analysis. 

Thermal and moisture effects on the strain energy release rates for interlaminar 
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fracture of unidirectional graphite/epoxy have been investigated by Russell and Street 
[13]. This investigation also included a study of the effects of shear loading through 
the use of various test configurations (Double Cantilever Beam, Cracked Lap Shear 
etc.). Initiation energies for delamination were found to increase as the proportion 
of shear loading increased and as the temperature was lowered, but no significant 
moisture influence was observed. The fracture resistance to crack extension was found 
to increase under tensile dominated loadings with both temperature and moisture 
content, but for high shear loading, the resistance was insensitive to the hygrothermal 
conditions. 

O’Brien, Raju and Garber have presented a CLT based analysis of mixed mode 
edge delamination specimens including hygrothermal effects [14]. They have used 
finite element modeling to determine the strain energy release rate components. Their 
results indicate total strain energy release rate increases of as much as 170% due to 
thermal effects for some T300/5208 graphite/epoxy laminates. However, a moisture 
content of 0.75% has been shown to totally alleviate this increase. According to 
. this analysis, in general, the consideration of thermal effects increases the energy 
release rate w’hereas moisture effect s have the opposite influence. These results have 
been confirmed using shear deformation models in the case of off-mid-plane edge 
delaminations [15]. It was found that the interlaminar stresses follow the same trend 
as the energy release rate, with increase due to thermal effects and alleviation due to 
hygroscopic effects. 

Aoki and Hondo calculated the strain energy release rate under thermal loading 
for mixed mode edge delamination. They used conventional finite element method 
[16] and a simplified method [16, 17] based on the classical lamination theory in com- 
bination with the J-integral for mechanical loading. Two types of axial constraint 
conditions were considered : (1) constant strain or fixed-grip and (2) constant load. 
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Numerical examples for cross-ply and angle-ply laminates showed that in angle-ply 
laminate, the energy release rate under free axial elongation increased constantly 
with delamination growth, while it remained constant under fixed-grip conditions. A 
higher order plate theory including transverse normal strain and thermal effects has 
been developed by Whitney [18] for the analysis of mid-plane edge delaminations. 
This approach provides the interlaminar stresses also, in addition to the strain en- 
ergy release rate. A [ 03 / 903 ], graphite/epoxy mode I specimen was analyzed and 
the maximum interlaminar normal stress was shown to increase by a factor of 2.7 
due to thermal effects, when compared with the pure mechanical strain reference 
configuration. 

From this summary it is found that there is a need for a unified approach that 
includes the analysis of free edge as well as local delamination and their interaction. In 
practical composite configuration free edge delamination does not occur in isolation, 
it is accompanied by other damage modes. Developing an analysis methodology that 
includes the interaction of delamination with other damage modes is essential for 
designing damage tolerant structures. 

The study of delamination consists of two main sections. These are the analysis of 
mid-plane edge delamination and local delamination in laminated composite plates. 

2.2 Mid-Plane Edge Delamination 

A mid-plane edge delamination specimen is shown in Fig. 2.2. A uniform axial strain 
e is applied in the x direction. From symmetry only one quarter of the specimen is 
considered. The sublaminate scheme and the choice of coordinate axes are illustrated 
in Fig. 2.3. 

Sublaminates 1 and 2 in Fig. 2.3 represent the uncracked and the cracked regions, 
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respectively. The analysis is based on the following displacement fields within each 
sublaminate 

« = xe + U(y) + z/3 x (y) 

v = V(y) 4 zBy(y) (2.1) 

w = W(y) 

where u,t>, and tr denote the displacements relative to the x,y, and c axes, respec- 
tive!}'. Shear deformation is recognized through the rotations 0 r and (3 y . In the 
present formulation thickness strain is neglected, and consequently inaccurate values 
of interlaminar peel stress, <r I2 , are expected. However, the peel stress can be modified 
by enforcing the free edge boundary condition associated with the transverse shear 
stress resultant. 

A generic sublaminate along with the applied forces and moments is shown in 
Fig. 2.4. The force and moment resultants are denoted by A’,, Q „ and M,. respectively. 
The constitutive relationships in terms of these force and moment resultants can be 
written as 

Nt = A,jCj + J S tk K k - ( i,j . k = 1,2,6) 

Mi = Bijtj + D lk K k - Mr (iJ,k = 1,2,6) (2.2) 

Qi ~ AjjEj (iij = 4,5) 

The subscripts x,y,z,yz,xz, and xy are replaced by the subscripts 1-6, respec- 
tively. The non-mechanical forces and moments resulting from hygrothermal effects 
are labeled with superscript nm for non-mechanical. They are defined as 

+ - 

(A7"",MD = j ' {ATa i + AH? i }Q ij (l,z)iz 
J ~2 


(2.3) 
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The thermal coefficient is denoted by aj in Eq. (2.3), while the swelling coefficient 
by fiy The are the plane stress sublaminate reduced stiffnesses [19]. The bars on 
Qj , fij and Qij indicate that theses quantities are to be obtained through appropriate 
coordinate transformations. The change in temperature between the ambient and the 
stress free temperature is denoted by AT. The percentage moisture weight gain is 
represented by AH. 

For a sublaminate of thickness h. the elastic stiffnesses and Z?,j in Eqs. 

(2.2) are defined as 

h 

{Aij.Bij.Djj) = [ +2 Q tj (l,z,z 2 )dz (2.4) 

J ~2 

The equilibrium equations can be written as follows 

■^W.v + *2* ~ hx — 0 

Ny.y + tiy — Uy = 0 

Qv.y + V 2 ~ Pi = 0 (2.5) 

Mx y.y ~ Qr + £(*21 + tjx) = 0 
My.y — Qy+ £(^2 y + *1 y) = 0 

where t 2 r,t 2 y ,P 2 and t\x,ti y ,Pi denote the interlaminar stress components at the 
sublaminate upper and lower surfaces, respectivelj'. These stress components appear 
in Fig. 2.4. Partial differentiation is denoted by a comma in Eqs. (2.5). Application 
of the boundary conditions and the governing equations to each of the sublaminates 
results in a system of differential equations which are solved to obtain the stresses 
and strain energy release rate. The boundary conditions to be prescribed at constant 
values of y, the sublaminate sections, are N x y or U,N V or V, Q y or W,M y or /3 y and 
M xy or 
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2.2.1 Uncracked Region (Sublaminate 1) 

From symmetry conditions at the sublaminate bottom surface, both w and the shear- 
ing stresses are zero. Since thickness strain is neglected, this leads to w being zero 
everywhere in this sublaminate. The equilibrium equations can be written as 

A T *» ,y = 0 

Ayj ,y = 0 

Q*.v-Pi=0 ( 2 - 6 ) 

,y ~ 0 

~ Qyi = 0 

where subscript 1 refers to sublaminate 1. From Eqs. (2.6) and the continuity of axial 
and in-plane shear stress resultants between sublaminates 1 and 2, we get 

N V1 = N xyi = 0 (2.7) 

By substituting from the constitutive relations into Eqs (2.6) and Eq. (2.7), and 
assuming an exponential form for the rotations @i y and (3 ix . we get the "following 
characteristic equation 

E a s 4 - E 2 s 2 4 E 0 = 0 (2.8) 

with 

Eo = A44^S5 — (j445) 2 

E 2 — O21 j4s5 4 ^32-4.44 — 0,22 — ^31^45 
E 4 — fl21^32 — ^22^31 

where 


(■^12^12 4 •'4-16^22 4 B 12 ) (j4i 2^13 + -Al6$23 4 B\$) 

[fl] 3x2 = (^22^12 4 B 2 e(22 4 D22) {B22U3 4 5 2 6$23 4 D 2 q) 

.(•526^12 + B^22 4 B 2 &) {B 2 6£l3 4 ^66^23 4 D$s) m 
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4 2 2 

■4 2 6 

-J 

_ -4 2 6 

-466, 



An A26 An B22 B26 

lsJix3 = ~ . _ _ 

.>*26 "66 J - 4 i 6 ^26 -£> 66 , 

Coefficient E 0 depends solely on the sublaminate axial stiffness, while E 4 is pre- 
dominantly influenced by the bending and coupling coefficients D,j and B,j . The 
numerical value of E 4 can be orders of magnitude smaller than E 2 and Eo ■ This results 
in the presence of a boundary zone in the response. For the material and laminate 
layups considered, the roots of this characteristic equation are real. Only the negative 
roots of Eq. (2.8) are considered as they give solutions decaying exponentially from 
the delamination tip. The solution can be written as 


where 


Ije-* v 


__ ^21 Sj ~ A 44 

^22-sj — -445 


(0 < y < b — a) 


0 = 1,2) 


Parameters Ij are arbitrary constants to be determined from the boundary condi- 
tions. By substituting Eq. (2.7) into the constitutive relations and using the assumed 
displacement fields, we obtain 


{*} + 


62 63 J 

(22 Z 23 I fill, 


( 2 . 10 ) 


where 


f5,j _ A 22 A 26 j Ny r- 
\ $2 j . -426 4.66 _ \ A*„ J 

Substitute from Eqs. (2.10) into the constitutive relations to get the resultant forces 
and moments in terms of the total extensional strain 

A*, -4n 4 -4 12^1 1 4 4i6^2i 

' M Vl > — Bn + -B 22 £n 4 B 26 C 21 {e} 

, M xy , , . 6 4 B 2 661 4 ^66^21 . 
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'A 2 sr + ^ 6 sr - K m i 

4 B 22 sr 4 BnS™ - M™ 4 [fi] 3x2 r lw,v 1 (2.11) 

. B 26 sr 4 BeeSr ~ M™ J ^ 

2.2.2 Cracked Region (Sublaminate 2) 

From the stress free boundary conditions at the face y = —a of sublaminate 2 and 
the equilibrium equations, we get 

■^1/2 = XV2 ~ Qv I = — o 

The equilibrium equations reduce to 

~ Qr 2 = 0 

Following a similar procedure as in sublaminate 1, the rotation can be written as 

(3 2x = Hi€ Xy (~a<y< 0) (2.12) 
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terms of fii< nlal strain 


where 


I v --- I . 

I M. m I 


{0 


4 


Ayy 4 Ai 2 <r > H + -^4.16^21 4 By 2 ip 3 y 

B 16 + ^26^11 4 B66<P2\ 4 D 2 ^<p 3 y 
A 12 Fr 4 A 16 Fr 4 B 12 F? m - A 

B 2(t F™ 4 BeeF™ 4 D 2 «Fr - M. 
Bye 4 j4j2^12 + Ay$<p 22 4 By 2 ^> 32 

Dee 4 B 2 6<Pl2 4 B 66^22 4 D 26V J 32 


Tnm 

a* 


nm 


AA,e Ay 


Fy' 

nm 

By ' 

f 2 

«— « 

1 

> 

11 

N X y 

.f 3 . 


.My . 


(2.13) 


The response associated with sublaminates 1 and 2 shown in Fig. 2.3 is coupled 
through the following conditions at y = 0. 


M Vl (0) = M w (0) = 0 
M iyi (0) = M xw (0) = 0 
/MO) = MO) 

The solution for both sublaminates i.e., the values of lj and Hy can be obtained 
by applying these conditions. The final expressions for the sublaminate rotations is 
given by Eqs. (2.9) and Eq. (2.12) where 

_ {-©i + (©4 + © 5 ^ 2 )( )} £ + ©s’M^) - ©2 4 ©4(f^) 

©3 — (©4 4 © 572 )^ 4 ©5»?1 

I 2 = --^{Aye 4 A 2 4 A 3 / 1 ) 

By = l/Ji 4 T] 2 h 

with 


Ai = By 2 4 B2261 ■+ B 26 ( 2 y 

a 2 = B 22 sr 4 b x st - Ar; m 


18 


As — —(^21 4 - Vl^22)si 
A4 = —(^21 4 ^2^22 ) s 2 

and 

©1 = i?26(£ll — ^ll) 4 ^66( ^21 ~ <p2l) ~ D26<P21 

©2 = B 26 (sr - Fr ) 4 Bee(Sr - F” m ) - DnFZ m 

©3 = — (^31 4 »h^ 32 )^l 

©4 = — (^31 4 T? 2 ^ 32) s 2 

©5 = ~(D 66 4 -626*1212 ■®66*r’22 4 D 2 $ ips 2 )\ 

The total strain energy release rate can be calculated by considering the work 
done by the external forces. This is given by 

2 jti r 

Gr = C, = -l < 2 - 14 ) 

where W, = work done by the external force in sublaminate i.L — laminate length, 
and a = crack length. 

The work done by the external forces is written in terms of the mechanical strain. 
£?, as 

W, = § [ e?N x ,dy (2.15) 

J v> 

with 

c” = £ - e”" 1 (2.16) 

The sublaminate free expansion strains, e" m , are calculated by setting the axial force 
resultant to zero, i.e. 

f N' i dy = 0 


(2.17) 
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nm 
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(AT+AH) AT=AH=0 
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N 


x 


0 



+ 


nm 

£ * 



C 


(a) 


(b) 


(c) 


Figure 2.5: Effective non-mechanical free expansion strain across the entire width of 
the laminate 


where N Xt is given in Eq. (2.11) and Eq. (2.13). The expression for each sublaminate 


is 

A C nrn _l A C nrn \rnm 

e nm _ ^12Q 3 "t~ A ie b 2 - A g 

1 An + -4i2^u + ^4j6^21 

£ „ m = A u Fr + A, e Fr + B 12 Fr - NT 

An + Anipn + A\^2\ + Bu'Pzi 

The total strain, £, is given by 


(2.18) 


e = e m 4 e nm 


(2.19) 


where t m is the effective mechanical strain and e nm is the effective free-expansional 
strain across the entire width of the laminate estimated by decomposing the non- 
mechanical problem in Fig. 2.5(a) into the superposition of the two cases shown in 
Fig. 2.5(b) and Fig. 2.5(c). In case (b) the laminate is subjected to a non-mechanical 
change ( AT + A H), while the strain is prescribed to be zero. In case (c) a unit axial 
strain is applied, while no hygrothermal change is considered. Knowing that no axial 
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force is applied in the main problem, i.e. Fig. 2.5(a), the sum of the axial forces in 
the two subproblems should be zero, hence 


N Xb 4 c nm N Xc = 0 


( 2 . 20 ) 


nm A T rfc 

~ ~ A T cc 


( 2 . 21 ) 


where N Xi and N Xc are the axial forces in case (b) and case (c), respectively. These 
axial forces are computed by substituting the expressions of N X) and N Xi from Eq. 
(2.11) and Eq. (2.13) into the relations 



r(b-a) 

r° 

11 

•o 

/ A*, dy 4 

f N x ,dy 


Jo . 

t J, 

f 

/■(fc-o) 

r° 

II 

/ A r x dy 4 j 

f N* s dy 

l 

Jo J 

Jr 


0,(A T+AH) 


The expressions for N Xb and N Xc are found to be 

N Xb = 4 A i6 Sr - ND(b-a) 

4 (A l7 Fr 4 A^F™ 4 B U F™ - N™ )o 

and 


( 2 . 22 ) 


(2.23) 


N Xc = ( Aii 4 A\ 7 in + -Aisfoi )(h — a) 

+ (^n + A^'Pw 4 ^16^21 + Bi2<P3i)a (2.24) 

The crack length and half of the total laminate width are denoted by a and b, respec- 
tively, as shown in Fig. 2.3. 

By combining the expressions of N Xi from Eqs. (2.11) and Eqs. (2.13) with Eqs. 
(2.14)-(2.24), the total energy release rate for the Mode I case can be written in the 
form 


Gj = (Gjl) + (Gil .) ■+ (Gin 4 Gir„ ) 


(2.25) 
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where 

Gil — (£ — C] m )[(- 4 ii 4 - 4 -^16^21)^ 4 4 AieS™ — N”™)) 

— ( £ ~ £2™ ) M11 4 Ai2<pn 4 - 4 . 16 SP 21 4 £12^31 ) £ 

4 (A 17 Fr 4 A ie Fr 4 BnF™ - A™)] 

GlLa — — ~ a ) {Mil 4 -^1261 4 - 416^21 ) [( 2 t — e" 7 ”) 

4 (Avsr+Ansr-jvr)]} 

dc 

— —a {(An 4 + ^le^i 4 -612^31 ) (2r — e 2 m ) 

. + (Ai 2 F? m 4 A 16 F 2 nm + - K m )} 

Gir — — (e — e”" 1 ) [(fin 4 171^12) I\S\t * l(i> o) 4 (fin 4 172^12) h s 2 e * 2<i> o) ] 

— (£■ — e”™ ) (Bi 6 4 j 4 j 2 'r > 12 4 v 4 i 6 i,P 22 4 B12P32) H^Xe Xa 

Gir 0 — — — {(^11 4 171^12) h [ e " l(i> a) — lj 4 (fin 4 172^12)^2 [f * j(fc °* — l] 

4 (i?i6 4 w4i2v?i2 4 Aie ( r > 22 4 -612^32) Hi — e *“) | 

and 

de_ _ ds™ _ ^12(5^ - F™) 4 A ie (S^ m - F 2 m ) - B 17 F™ 

da da N Xc 

Al6 l>il2(^ll ~ <pll) + j4ie( ^21 — <^2l) — Bi2<p3l ] ln 0/ , N 

\T2 (l.£0) 

J x< 

The resulting expression for the total energy release rate Gj is composed of three 
terms. The first term, denoted by Gjl is independent of the delamination length while 
the second, Gu a , is a linear function of the crack length. The third term denoted by 
(Gir 4 GiR a ) , is an exponentially decaying function of the delamination length. 

In computing the non-mechanical strains, the laminate is assumed to be held at 
the prescribed temperature and moisture levels. This is followed by testing under 
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fixed-grip condition, i.e., the constant strain measured in the lab is the mechanical 
strain e m . In Refs. [18] and [20], Whitnej' considered the strain measured in the tests 
to be the total strain, i.e. e = c m + s nm = constant . The difference between the 
two interpretations is detected by the terms Gn a and Gjr b in Eq. (2.25). These two 
terms are neglected in Refs. [18] and [20] since the total strain e is assumed to be 
constant. 

As mentioned previously, neglecting the thickness strain leads to inaccurate esti- 
mates for the peel stress. The peel stress is given by 

V = Q Vl ,y = -(A 4A + A A sTij)IjSj€~*’ v (j = 1.2) (2.27) 

The equilibrium of transverse force requires that 

/ pdy = 0 (2.28) 

Jo 

or from the equilibrium equations (2.6) 


Q 


Vi 




- Q 




= 0 


While for all practical purpose the resultant shear stress Q y vanishes due to 

the free edge, the resultant shear stress at the delamination front Q y =£ 0. That 
is in order for the peel stress to satisfy transverse force equilibrium, the shear force 
boundary condition at the sublaminate end should be enforced. This is done by 
adding to the peel stress distribution an appropriate boundary function expressed in 
terms of the characteristics roots as 


a\e~ tlV -j- ait~‘ 3V 


The coefficients aj and a 2 are obtained by enforcing equilibrium of transverse force 
given in Eq. (2.28) and moment given by 



-f M Vi 


*»i =(»-•) 


= 0 
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Figure 2.6: Local Delamination Specimen Cross Section 
The corrected peel stress distribution is 


P = 


S1S2 


-M, 


s ] — s 2 Vi l »i=<»-») 


[ 53 e-*»v_ 52e -.,y] 


2.3 Local Delamination 


A longitudinal section illustrating the geometry of a generic configuration is shown 
in Fig. 2.6. The central region is assumed to be made of 90° plies with an isolated 
transverse crack in the middle. Delaminations are assumed to grow from both ends 
of the transverse crack, and towards both specimen ends as shown. From symmetry 
considerations, only one quarter of the configuration is modeled. The modeled portion 
of length L is divided into four sublaminates as shown in Fig. 2.7. The crack length 
is denoted by a. The top surface (sublaminates 1 and 4) is stress free. In order to 
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free free 



Plane w=0 w=0 


Figure 2.7: Sublaminate Scheme for Local Delamination 

simplify the analysis, the thickness strain e z is neglected. The consequence of this, 
combined with the fact that the transverse displacement w is zero along the center line, 
is that w is zero in sublaminates 1, 2, and 3. Also, this approximation does not allow 
for the enforcement of boundary conditions on the shear stress resultants, leading 
to incorrect estimates of the interlaminar normal stresses. The interlaminar shear 
stress estimates, however, are reliable [6]. These assumptions lead to considerable 
simplifications in the analysis. In spite of the simplifications, reliable energy release 
rate components can be estimated based on the interlaminar shear stress distributions 

[7]. 

A generic sublaminate is shown in Fig. 2.8 along with the notations- and sign 
conventions. The peel and interlaminar shear stresses are denoted by P and T, 
respectively, with t and b subscripts for the top and bottom surfaces, respectively. 
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Figure 2.8: Generic Sublaminate for Local Delaxnination 

The axial stress resultant, shear stress resultant, and bending moment resultant are 
denoted by A T ,Q, and Jlf, respectively. The governing equations correspond to the 
one-dimensional form of Eqs. (2.1 - 2.5). These are summarized in the following for 
convenience. 

The x and z displacements within the sublaminate are assumed to be of the form 

u(x,z) = U{x) + z/3{x) 
w(x,z) = H T (x) 

Here, U represents the axial mid-plane stretching and W is the transverse displace- 
ment. The shear deformation is recognized through the rotation, /3. These displace- 
ments are the total quantities and include the hygrothermal effects. The origin of 
the coordinate axes for the sublaminates is taken at the delamination tip as shown 
in Fig. 2.9 . The equilibrium equations take the form 

+7) - Tb = o 


Qy . r +Pt — Pb — 0 


(2.29) 
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Figure 2.9: Sublaminate Forces and Coordinate Systems 
M ix -Q + ^(T, + T b ) = 0 

where h is the thickness of the sublaminate. The constitutive relationships in terms 
of the force and moment resultants are 

N = A U U, X 

Q = An(/3+W,?) 

M = B n V, x 

The boundary variables to be prescribed at the sublaminate edges are 

N or U 
M or 0 
Q or W 
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Additionally, at the interfaces between sublaminates, reciprocal traction, and dis- 
placement matching boundary conditions have to be specified. The stress resultants 
in these equations include the equivalent hygrothermal loads also. 

The solutions in sublaminates 1 and 2. are coupled by the reciprocal interlaminar 
stresses denoted 7\ and P\ and by displacement continuity at the common interface. 
Assuming exponential solutions for the axial force and bending moment resultants 
leads to an eigenvalue problem involving the exponential parameters. The character- 
istic equation is of the form 

-j- B%s 2 -f .63] = 0 

where s is the eigenvalue parameter, and the B coefficients are given by 

Bi - ( * + 1 + + ^2 \ AilJJ Ai(2) 

1^11(2) 4D n (i) 4 Z>h(2) / A 55( i) A 55(2) 



For the problem under consideration, all the square roots in this expression lead to 
real quantities and thus the eigenvalues are real. Since the eigenvalues involve only 
the stiffness parameters, they are not affected by the inclusion of hygrothermal effects. 
Further, due to the fact that B^ has D terms in the numerator, it is much smaller 
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than £3. This leads to the boundary layer nature of the solution. Since the response 
(axial forces, moments) has finite values at large distances from the origin, namely, at 
the ends of the specimen, only the exponentially decaying and constant solutions are 
used. Using subscripts to denote the sublaminate of validity, the following boundary 
conditions from the ends of the modeled region are enforced. 

JV a (0) = 0 


QM = 0 

fata) = 0 

A T i -f A r 2 = Applied Load 

The conditions on N apply only to the mechanical quantities. Further, the fol- 
lowing displacement matching conditions are applied. 



Vi(0) = U 4 ( 0) 
V 2 (0) = U 3 (0) 
01(0) = 04(0) 


It should be noted that a 0 2 and 0$ matching condition cannot be applied at this 
level of modeling since it would amount to specifying both W and Q. To eliminate 
rigid body displacements, U\ is set to zero at the left end. The following solutions 
can then be obtained for the stress resultants in sublaminates 1 and 2 


Ny = oie* 1 * + a2e* JI -1- cAji(i) — A r ” m 
A 2 = — aje * 12 — U 2 e * SI + £^11(2) — 
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M, = c,fcie* ,x + a 7 k 2 e ,3t - M™ 

M 2 = a jfcae* 1 * + o 2 fc 4 e‘’ x - M™ 

Here ki is defined as 

*1.2 

fc, = 7 - 2 - 1 - 

•"til 1 ) _ .2 
■^ll(l) 1 

The parameter k 2 is defined in a similar manner using the eigenvalue, s 2 . The re- 
maining parameters, k 2 and fc 4 , are similar to fcj and k 2 but based on sublaminate 2 
properties. The nominal strain, £, is defined as 



where P is the applied uniform axial force and b is the specimen width. The a' s can 
be derived from the boundary conditions as follows 



— $3 — $1 + ($4 ~ $2)° 
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The interlaminar shear and pet*] stresses between sublaminates 1 and 2 can be ob- 
tained using the equilibrium equations (2.29) as 


Ti = fl.jSi f * ,T -j- O-zSzt'** 

P\ = (ki + — + (A'2 + — )fl2S2f* jX 

As mentioned previously, this peel stress estimate is not accurate because of the 
inability to apply boundary conditions on shear. Recognizing the fact that there 
are no applied shear forces, it can be concluded that the peel stress distribution 
should be self equilibrating. This assumption can be satisfied by including additional 
exponential terms in the above peel stress expression and determining these additional 
terms by setting the net force and moment due to the peel stress to zero as shown 
in section 2.2. The peel stress estimated through this correction process is referred 
to as the modified peel stress. Proceeding on to sublaminates 3 and 4, the following 
solutions can be written. 

A 3 = 0 


where 


and 


Mz = sinh(w 3 x) -f ipi cosh(u> 3 x) 


lf>2 — a 1^3 "t" a 2^4 


COth(u>3<l) 



^55(2) 

■^ 11 ( 2 ) 


A T 4 = 


P 

2b 


M 4 = Qi&i -f a 2 fc 2 
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The total energy release rate Gj is calculated using Gt = dW e /da where W’ £ is the 
work done per unit width bj r the external (constant.) loads on the specimen displace- 
ments. For the case where hygrothermal effects are included, there are additional 
terms due to the work done by TV,"™ . In reality, these A 7 "™ quantities are not applied 
loads but correspond to residual stresses. Thus, the additional terms are due to the 
work done by the applied mechanical strains on these residual stresses. The total 


energy release including hygrothermal effects is given bj r 


+ IL N nm ( 1 

2 b 2 \-A ll{i) + 


•^■ 11 ( 1 ) -^ 11 ( 1 ) + ^ 11(21 


+ h ~ h 


— I 3 + /; 


ll(i) -r ^ii(2) 


where the I factors are 


h = x- 


/ 2 = X 


e 2 e 3 - e l e 4 (\ - i - e -‘* {L - a) 


(0 3 + e 4 a)e- ,i{L ~ 0) - (0, + 0 2 a)c- 3{L - a] 


■^ll(l) + A.1U2) >lll(l) 

Parameter 1 3 is the same as but with the ratio -An(i)Mii( 2 ) instead of unity in 
Eq. (2.31). Using the virtual crack closure technique [21], from the relative displace- 
ments in the cracked portion and the interlaminar stresses ahead of the crack tip, the 
mode I and mode II energy release rate contributions can be obtained. The mode III 
energy release rate is zero from the assumption of plane strain. The mode II energy 
release rate is given by 


Gu = lim f Ti(x - 6)Au(x)d2 

6—*0 ZO JO 


where 6 is the virtual crack step size and Au is the differential axial displacement 
across the crack surface. This cdculation can be simplified using only the linear 
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pari of the differential displacement [7]. In a similar fashion, the mode I energy 
release rate can be obtained based on the normal stress ( P ) and the differential 
tt> displacements near the crack front. Since the unmodified peel stress estimate is 
inaccurate, an alternate approach was used to estimate Gj, the mode I energy release 
rate. The total energy release rate for this problem is made up entirely of Gj and Gjj 
( Giii — 0). From an estimate of Gj and G/;, an estimate for Gj can be obtained 
simply as 

Gj = Gj — Gjj 

The critical load for a given specimen can then be evaluated based on an appropriate 
fracture law. This is illustrated in the next chapter. 
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CHAPTER III 

APPLICATIONS OF DELAMINATION MODELS 


3.1 Mode I Edge Delamination 

The anafytical model is applied to the mid-plane edge delamination specimen shown 
in Fig. 2.2. The material considered is T300/5208 graphite epoxy. Its properties are 
listed in Table 3.1. 

The difference between the ambient and cure temperature, AT, is —156°0. The 
moisture level was allowed to vary from 0 to 1.2 percent of the laminate weight, which 
reflects feasible conditions. Laminates of the class \9j— 02/0/9O 2 ] 4 and [0 3 /90 3 ] 3 have 
been analyzed. 

Normalized values of strain energy release rate are shown in Figs. 3. 1-3.6, where 
the labels M , M 4- T, and M -{■ T 4 H stand for mechanical, mechanical and ther- 
mal, and mechanical, thermal and moisture, respectively. The strain energy release 


Table 3.1: ED Specimen Geometry and Material Properties 


E n = 128 GPa 

Thermal Coefficients : ctj = — 0.41p£/°C 

E 22 = 8.47 GPa 

Q 2 = 26.8fie/°C 

G \ 2 — 5.73 GPa 

Swelling Coefficients : & = 0 

G 3 i 3.2 s GPa 

fit = 5560 fi£/%W 

(? 2 3 — 3.2 s GPa 

Width = 2 b = 38.4 mm 

i/ 12 = 0.292 

Ply Thickness = 0.14 mm 
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rate parameter in the figures is defined as the total energy release rate divided by 
E 22 h(e m )\ 

The strain energy release rate in Figs. 3. 1-3.3 is zero at a = 0. Residual thermal 
stresses results in an increase of 275%. 40% and 280%. of the energy release rate for the 
[15/ — 15 2 /15/90 2 ],, (60/ — 60 2 /60/90 2 ], and [03/90 3 ], laminates, respectively. Residual 
moisture alleviates this effect as illustrated in Figs. 3.4-3.6. The specific moisture 
content for total alleviation from the thermal effect is equal to 0.763%. irrespective of 
the layup. 

The peel stress distribution, a :: , appears in Figs. 3. 7-3. 9. The stress parameter in 
these figures is defined as the interlaminar stress divided by E 2 i£ m • The inaccurate 
peel stress distribution given in Eq. (2.27) is plotted for the case where mechanical 
loading only is considered. The corrected peel stress distribution is self-equilibrating 
and yields a tensile peel stress at the delamination front. 

The magnitude of the peel stresses shows a strong dependency on the thermal 
and moisture conditions. The stress increases with thermal effect as compared to 
pure mechanical loading.. The addition of moisture alleviates the thermal effect. 
Moreover, the distance at which the peel stress reverses its sign is not affected by the 
residua] thermal and moisture strains. It is worth noting that at the specific moisture 
percent (0.763%) producing complete alleviation of the total energy release rate from 
the thermal effect, the interlaminar peel stress distribution is identical to the case 
where only mechanical loading is considered. This is shown in Figs. 3. 7-3.9. This 
finding establishes a similarity in behavior between the energy release rate and the 
interlaminar stresses. 

The analytical model presented herein was applied to the laminates presented in 
Ref. [18]. The Graphite/Epoxy l ami na properties from Ref. [18] are listed in Table 
3.2. S imil ar values of strain energy release rate Gi were calculated for the wide range 
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Table 3.2: ED Specimen Geometry and Material Properties, Ref. [18] 


E 11 /E 22 — 14 

Ply Thickness = 0.1267 mm 

E 33 /E 22 = 1 


G 12 /E 22 — 0.533 

Thermal Coefficients : Qi = —0.9 ns/°C 

G 23 /E 22 = 0.323 

q 2 = ct 3 = 23.0 fi£/°C 

V \2 — 0-3 


V 23 = 0.55 

width = 2 b = 38.0mm 


of a/h where the Gj remains constant as shown in Figs. 3.10 and 3.11. Negligible 
change in the Gi value with decreasingly small values of a/h were obtained. This 
is in contrast with the increase in Gi at small values of a/h reported in Ref. [18]. 
Although thickness strain is neglected in Eqs. (2.1), the peel stress distribution has 
been estimated through a modification as described previously, which simplifies con- 
siderably the computational effort . A comparison of the peel stress distribution with 
Ref. [18] is shown in Figs. 3.12 and 3.13. 

The peel stress intensity at the delamination front in the [30/ — 30 2 /30/90 2 ] 4 is 
higher than the [0 3 /90 3 ], laminate. This is due to the difference in poisson’s ratio 
between the core plies made of 90° plies and the outer plies. The poisson’s ratio 
mismatch is larger for the case of [30/ — 30 2 /30/90 2 ], compared to the [0 3 /90s], layup. 
The interlaminar peel stress distribution predicted by the present approach is in good 
agreement with the distribution of Ref. [18] for the case of a [0 3 /90 3 ],laminate. This is 
in contrast with the case of a[30/ — 30 2 /30 /90 2 ], where the maximum stress intensity 
as well as the distribution differ from the predictions of Ref. [18]. This difference may 
be due to the transverse normal strain influence on the analysis of these laminates. 
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3.2 Edge and Local Delamination 

The delamination models have been used to study the behavior of [±25/90 n ], T300/934 
Graphite Epoxy specimen for n values of 0.5, 1, 2, 3, 4, 6, and 8. These correspond 
to the specimen tested by Crossman and Wang [8]. The specimen width and length 
were fixed at 0.025m and 0.15m, respectively, as in the tests. In computing the non- 
mechanical strains, the laminate is assumed to be held at the prescribed temperature 
and moisture levels. In predicting critical strains, the difference between test and 
stress free temperatures is assumed to be — 155°C' and specimen is assumed to be 
dry. It is assumed that local delamination occurs under fixed load conditions whereas 
edge delamination occurs under fixed grip conditions. This difference is a consequence 
of the modeling approaches used in the analyses. The applied uniform load was 100 
MPa axial stress for the local delamination analysis and 0.5% strain for the edge de- 
lamination analysis. The solutions were generated using simple computer programs 
based on the closed form expressions for the interlaminar stresses and energy release 
rates. 

3.2.1 Local Delamination 

An example of the total energy release rate variation associated with local delamina- 
tion (neglecting hygrothermal effects) with the crack length is presented in Fig. 3.14. 
The asymptotic value of Gj is denoted by Gto in the figure. It can be observed that 
after a certain crack length, the Gt is independent of the crack length. On the basis of 
curves like the one shown in Fig. 3.14, the crack length was fixed at 10 ply thicknesses 
for the remainder of the studies. Typical interlaminar shear stress profiles including 
the hygrothermal effect are presented in Fig. 3.15. The corresponding total strain 
energy release rates appear in Fig. 3.16. The inclusion of thermal effects increases 
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the stress and the energy release rate while the inclusion of moisture effects has the 
opposite effect. In fact a moisture level of about 0.75% almost exactly negates the 
thermal effects. After some initial dependence on crack length, the mode mix tended 
to stabilize to a constant value. Using the model developed here, the asymptotic 
mode II component of the local delamination energy release rate was found to be 
approximately 30 percent for all n values. In the case of off-mid-plane edge delami- 
nation, the mode II contribution was less than 10 percent for the n = 0.5 specimen 
and progressively less for the thicker specimen. 

3.2.2 Edge Delamination 

As in the case of local delaminations, the interlaminar stress increases with thermal 
effects and the addition of moisture alleviates this as shown in Fig. 3.17 for the case of 
mid-plane edge delamination. A moisture level of about 0.75%. produces a modified 
peel stress distribution that is indistinguishable from the case of mechanical loading 
alone. Moreover, the distance at which the modified peel stress reverses its sign is 
not affected by the residual hygrothermal strains. The hygrothermal influence on 
mid-plane delamination strain energy release rate is illustrated in Fig. 3.18 where the 
strain energy release rate is plotted versus moisture content for a [±25/902]* laminate. 
The strain energy release rate follows the trend of increasing with residual thermal 
stress as in the case of peel stress. Further, residual moisture alleviates the thermal 
effects and a moisture level of about 0.75% results in a total alleviation of thermal 
effects. Similar behavior is observed in the case of off- mid-plane edge delamination. 

3.2.3 Failure Loads and Modes 

In order to evaluate the critical loads for local delamination, an appropriate mixed 
mode fracture law has to be applied, based on the calculated energy release compo- 
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nents. The following mixed mode criterion [22] has been fitted to the test data of 
Ref. [23] to calculate the mixed mode Gt c which is then used in the Griffith criterion 
Gt = Gj c to obtain the critical delamination growth stress <r c and strain c c values. 

Gtc = CGic + (1 — iX'Gjic 

Here £ is the mode I fraction ( Gi/Gt ) and Gj c and Gjj c are the critical strain energy 
release rates for the limiting cases of pure mode I and pure mode II, respectively. 
The exponential parameter tj is a material constant and for the T300/934 system, its 
value is approximately 0.9. In the case of mid-plane delamination, since only mode 
I is present, Gt c was taken as G; c ( 125 J/m 2 ). Based on the mixed mode criterion, 
Gjc w as about 400 J/m 2 for the local delamination case (£ = 0.7). The failure loads 
for edge delamination at the -25/90 interface have also been calculated using a Gt c 
value of 150 J/m 2 . This Gr f value is different from the value used for mid-plane 
delamination due to the limited (less than 10 percent) presence of mode II. 

In order to consider a worst case situation, thermal stresses were included and the 
moisture level was set at zero. Though the thermal stresses had a significant effect 
on the calculated peak stresses, the effect on the energy release rate was not signifi- 
cant except in the case of mid-plane edge delamination for the material system and 
layup considered. The critical strains are plotted against n, the number of 90° plies 
in Fig. 3.19. The experimental results of Ref. [8] are also presented in the figure for 
comparison. The results of the model developed in this paper are represented by the 
solid and dotted lines while the experimental results are shown as filled squares. The 
CLT based model of Ref. [9] agrees well with the shear deformation model in terms 
of the total energy release rate. However, the CLT based model does not provide in- 
formation on the mode split and thus, the value of G c (fc Gj c ) used can lead to bias in 
the critical strain estimates. In the experiments, the local delamination phenomenon 
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was observed as the predominant failure mode only for the n = 4,6 and 8 specimens. 
The shear deformation model presented in this paper provides good agreement with 
the experimental data in this range. For n < 4, edge delamination either in the 
mid-plane or in the 25/90 interface was observed in the tests, in agreement with the 
edge delamination models. Further, the relative closeness of the calculated critical 
strains from the mid-plane and off-mid-plane edge delamination models implies that, 
in practice, one could have interaction between these two modes. In such cases, one 
can expect the delamination to wander around the mid-plane and the 25/90 inter- 
faces. This is especially so in the case v = 0.5 where mid-plane delamination is not 
actually between two distinct layers but in the middle of a single layer. Experimental 
observations [8] are in agreement with this expectation. Thus, it can be seen that 
the shear deformation models reproduce the observed behavior with reasonable ac- 
curacy and can be used to estimate critical loads for a range specimen thicknesses 
incorporating various delamination modes. 

3.3 Conclusions 

Shear deformation models including hygrothermal effects have been developed to 
analyze local delaminations growing from transverse cracks in 90° plies and edge 
delaminations located around the mid-plane of symmetric laminates. The models 
have been combined into a unified delamination analysis code in order to predict 
damage modes and loads in laminated composites. The analytical results of the 
shear deformation models agree reasonably with critical strain experimental data 
from [±25/90 n ], T300/934 graphite epoxy laminates in the range of n from 0.5 to 8. 
Residual thermal and moisture stresses are found to have only minor effects on the 
critical strains except in the case of mid-plane edge delamination for the geometry 
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and material considered. The same failure modes as in the tests are reproduced in 
the analysis. The integrated delamination code is expected to be of nse in design 
evaluation applications. 
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Figure 3.2: Mode I Strain Energy Release Rate in a [60/ - 60 2 /60/90 2 ], Laminate 











Graphite/Epoxy 

| 

[60/-60 2 /60/90 2 ] s 

l 

e m *=0.00254 

1 

a/h=10 

! 

i 


1 

I 0.763% 


2 0 . 


Figure 3.5: Influence of Residual Thermal and Moisture Stresses on Mode I Strain 
Energy Release Rate in a [60/ - 60 2 /60/90 2 ], Laminate 
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Figure 3.7: Peel Stress Distribution ahead of the Crack in a [15/ - 15 2 /15/90 2 ] 
Laminate 
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Figure 3.8: Peel Stress Distribution ahead of the Crack in a [60/ — 6 O 2 /60 /90 2 ] 
Laminate 
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Figure 3.9: Pee] Stress Distribution ahead of the Crack in a [0 3 /90 3 ] a Laminate 
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Figure 3.10: Mode I Strain Energy Release Rate in a [45/ — 452/45/902]* Laminate 
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Figure 3.13: Peel Stress Distribution ahead of the Crack in a [30/ — 30 2 /30/90 2 ] 
Laminate 
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Figure 3.16: Total Energy Release Rate for [±25/90„], Graphite/Epoxy Specimen 
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CHAPTER IV 

THEORY OF ANISOTROPIC THIN-WALLED BEAMS 


A variationallv and.asymptotically consistent theory is developed in order to derive 
the governing equations of anisotropic thin-walled beams with closed cross sections. 
The theory is based on an asymptotical analysis of two-dimensional shell theory. 
Closed-form expressions for the beam stiffness coefficients, stress and displacement 
fields are provided. The influence of material anisotropy on the displacement field 
is identified. A comparison of results obtained by other analytical developments is 
performed. 

A review of previous work is presented first, this is followed by a detailed develop- 
ment of the theory. Finally an analytical comparison of the displacement field with 
previously developed theories is provided. 

4.1 Review of Previous Work 

Elastically tailored composite designs are being used to achieve favorable deformation 
modes under a given loading environment. Coupling between deformation modes 
such as extension-twist or bending-twist is created by an appropriate selection of fiber 
orientation, stacking sequence and materials. The fundamental mechanism producing 
elastic t ail oring in composite beams is a result of their anisotropy. Several theories 
have been developed for the analysis of thin-walled anisotropic beams. An extensive 
review is provided in Ref. [26]. A number of issues relevant to the research undertaken 
in this thesis is highlited in the following. 
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A basic element in the analytical modeling development is the derivation of 
the effective stiffness coefficients and governing equations which allows the three- 
dimensional (3D) state of stress to be recovered from a one-dimensional (ID) beam 
formulation. For isotropic or orthotropic materials this is a classical problem, which 
is considered in a number of text books such as Refs. [52]- [59]. 

For generally anisotropic materials, a description of the major approaches is pro- 
vided in Refs. [24]- [49] . A number of ID theories have been developed in Refs. [27], 
[28j. [30j, [42], [43], and [46]. A discussion of the displacement provided in these works 
is presented in the analytical comparison section of this chapter. 

Missing from the review of Ref. [26] and all other current publications is the work 
of Reissner and Tsai in Ref. [27]. It presents an exact solution to the governing 
equilibrium, compatibility and constitutive relationships of shell theory. Closed as 
well as open cross-sections were considered. However, the authors left to the reader 
the derivation of the explicit expressions for the stiffness coefficients. This may be 
the reason for their work to have been overlooked. These expressions are important 
in identifying the parameters controlling the behavior and in performing parametric 
design studies. Furthermore, the explicit form of the displacement field helps evaluate 
and understand predictions of other analytical and numerical models. 

A number of assumptions were adopted in Reissner and Tsai’s development re- 
garding material properties such as neglecting the coupling between in-plane strains 
and curvature which can be significant in anisotropic materials. It is important to 
assess the influence of these assumptions on the accuracy. This has been done in the 
present work by using an asymptotical expansion of the shell energy. 

Mansfield and Sobey [28] and Libove [29] obtained the beam flexibilities relating 
the stretching, twisting and bending deformations to the applied axial load, torsional 
and bending moments for a special origin and axes orientation. Their analyses are 
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similar. Although they did not refer to the work of Reissner and Tsai [27], surprisingly 
when their analyses is applied to the special case outlined in Ref. [27], their stiffnesses 
coincide. However, one has to carry out details to show this fact. They adopted the 
assumptions of a negligible hoop stress resultant N tt and a membrane state in the 
thin-walled beam section. The special case in Ref. [27] refers to the one where classical 
assumption of neglecting shear, hoop stress and constant shear flow is adopted. 

A pertinent element in the analytical modeling development is the inclusion of 
section warping. The major difference among the various theories lies in the method- 
ology used to eliminate warping and consequently obtain a one-dimensional theory. 
The work of Refs. [30], [41], [42], [43], [44], [45], and [46] use the displacement field 
of thin-walled isotropic beams with shear deformation as the basis of their analytical 
development. In Refs. [42] and [46] the torsional rigidity is derived in terms of Clas- 
sical Lamination Theory in what the author described as a “practical manner’'. In 
Refs. [43] and [44] a shear correction fact or has been introduced in order to reduce the 
overestimated bending stiffness. This factor was derived for the case of pure torsion-; 
by using the virtual work method and enforcing compatibility. While this approach 
shows an improvement in predictions, it is problem dependent. Another modification 
was proposed in the finite element formulation of Ref. [38]. This formulation aims at 
minimizing the error associated with the neglect of bending-related warping in the 
theory of Ref. [30]. This modification was based on shear stiffness correction factors 
determined by numerical comparison of results with an MSC/NASTRAN solution of 
cantilevered beam configurations loaded transversely at the free end. 

This summary points to the necessity of addressing three fundamental issues. 
The first, is the effect of the material’s anisotropy on the displacement field and how 
to include its contribution in a consistent manner. No rigorous proof is provided 
to validate the assumed displacement fields in Refs. [30], [42], [43], [44], [45], and 
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[46] for beams made of anisotropic material as indicated b} f the various correction 
factors introduced. The second, is the significance of the shear deformation relative 
to the other contributions such as section related warping. The last is the accuracy 
of the membrane stress state assumption in thin-walled anisotropic beam sections. 
The present work addresses these issues by using an asymptotical expansion of the 
2D shell energy to derive the ID beam displacement field. As a consequence, the 
material’s anisotropy is accounted for in a consistent manner and the deformation 
modes that have a lead contribution to the energy emerge naturally. 

4.2 Coordinate Systems 

Consider the slender thin-walled elastic cylindrical shell shown in Fig. 4.1. The length 
of the shell is denoted by X, its thickness by h, the radius of curvature of the middle 
surface by R and the maximum cross sectional dimension by d. It is assumed that 

d << L h << d h << R (4.1) 

The shell is loaded by externa] forces applied to the lateral surfaces and at the 
ends. It is assumed that the variation of the external forces and material properties 
over distances of order d in the axial direction and over distances of order h in the 
circumferential direction, is small. The material is anisotropic and its properties can 
vary circumferentially and in the normal direction to the middle surface as well. 

It is convenient to consider two coordinate systems for the description of the stat e 
of stress in thin-walled beams. The first one is the Cartesian system x,y and z shown 
in Fig. 4.1. The axial coordinate is x while y and z are associated with the beam 
cross section. The second coordinate system, is the curvilinear system x.s and £ 
shown in Fig. 4.2. The circumferential coordinate s is measured along the tangent 
to the middle surface in a counter-clockwise direction whereas £ is measured along 


h(s) 

d 


Figure 4.1: Cartesian Coordinate System 
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the outward normal to the middle surface. A number of relationships have a simpler 
form when expressed in terms of curvilinear coordinates. A relationship between the 
two coordinate systems can be established as follows. 

Define the position vector t of the shell middle surface as 

r = xi x + y(s)T y + z(s)u 

where i x , i y , T : are unit vectors associated with the cartesian coordinate system x. y 
and z. Equations y = y(s) and r = z(s) define the closed contour T in the y, z plane. 
The normal vector to the middle surface n has two nonzero components 

n = n v (s)r v + n.(s)u (4.2) 

The position vector R of an arbitrary material point can be written in the form 

R = r -f (ri (4.3) 



(4.4) 
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An asymptotical analysis is used to model the slender thin-walled shell as a beam 
with effective stiffnesses. The method follows an iterative process. The displacement 
function corresponding to the zeroth-order approximation is obtained first bj r keeping 
the leading order terms in the energy functional. A set of successive corrections is 
added and the associated energy functional is determined. The process is terminated 
when the new cycle does not generate any additional terms of the same order in the 
energy functional. 

4.3 Shell Energy Functional 

Consider in a 3 D space the prismatic shell in Fig. 4.2. A curvilinear frame x. s, and 
£ is associated with the undeformed shell configuration. Values 1, 2 and 3 denoting 
z, s, and £, respectively are assigned to the curvilinear frame. Throughout this 
study. Latin superscripts (or subscripts) run from 1 to 3, •while Greek superscripts 
(or subscripts) run from 1 to 2. unless otherwise stated. 

The strain energy density of a 3 D elastic body is a quadratic form of the strains 

V = ±E ijk, e, j£kl 

The material properties are expressed by the Bookean tensor E'^ kl . Following the 
classical shell formulation of [60], [61], and [62] the through-the-thickness stress com- 
ponents <r t3 are considerably smaller than the remaining components <? a0 . Therefore 
we can set 

cr i3 = 0 (4.5) 

so that the strains can be written as 


Caff — Hal 3 CPa0 


(4.6) 
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where 7 ^ and p a p represent the in-plane strain components and the change in the 
shell middle surface curvatures, respectively. For a cylindrical shell these are related 
to the displacement variables by 


7n = 


dv\ 

dx 


dvi dt > 2 

2711 = a7 + 17 

8v 2 V 

7:1 = a 7 + 1 

d 7 v 

Pn = 


dx 7 


d 2 X' t 1 
Pu = dlfrr + 4R 


fdvi di> 2 \ 
v ds dx) 


d 2 v d / 1>2 \ 

= I? ~ a* \r) 


P22 


(4.7) 


ds 7 8 s\R 

where Vj, t > 2 and v represent the middle surface displacements in the axial, tangential 
and normal directions, respectively as shown in Fig. 4.2. These are related to the 
displacement components in Cartesian coordinates by 


t'l = Vi 


dy dz 

v 2 = U 2 — V 3 — 

ds ds 

dz dy 

v = u 2 - Ha- 
as as 


(4.8) 


where Uj, v 2 and V 3 denote the displacements along the x, y and z coordinates, 
respectively. 

The energy density of the 2D elastic body is obtained in terms of -y a p and p a p by 
the following procedure. 

The 3 D energy is first minimized with respect to e » 3 . This is equivalent to satis- 
fying Eq. (4.5). The result is 

V = min C7 = ]rD afiyt e a pC^ 

f»3 2 * 


(4.9) 
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where D a ^ represent the component of the 2D Young’s moouhi ... The expressions 
for D al3 ' l( ' are given in terms of E by 


r»a/333 pyb 33 

D a0jt = 


£3333 


(4.10) 


where 


qo0h £>q0h2 _ 


£c»033£*i333 

£3333 


and H u \ are components of the inverse of the 2D matrix |j£ l ' 23Xa — 1|. The 

expression for in terms of familiar Classical Lamination Theory (CLT) param- 

eters is provided in Eqs. (4.43) and (4.44). 

The strain £ a p from Eq. (4.6) is substituted into Eq. (4.9). After integration of 
the result over the thickness ( one obtains the energy of the shell $ per unit middle 
surface area 

h 3 




2 / 100 -)/' 


a 0i6 


(4-11) 


where 


1 

h 


£10016 __ J-)O0i6 ^ 


Cf' t = f- < 


a3-)6 { 


lO0lf 


h 3 

= ~ < D aih *C > 


and a function of (, say a(£), between pointed brackets is defined as an integral 
through the thickness, viz., 

r + M *)/2 


a{()d( 

■ M «)/2 


(4.12) 


The first term in Eq. (4.11) represents the in-plane contribution, the second the 
coupling between in-plane and bending, and the third the bending contribution to 
the shell energy. 
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For an applied external loading P t , the displacement field u,- determining the 
deformed state are the stationary points of the energy functional 

I = J $dxds — J PiUidxds . (4.13) 

4.4 Asymptotical Analysis of the Shell Energy Functional 


4.4.1 Zeroth-Order Approximation 


Let A and E be the order of displacements and stiffness coefficients respec- 

tively. Assume that the order of the external forces is 



(4.14) 


This assumption is shown later to be consistent with the equilibrium equations. 
An alternative would be to assume the order of the external force as some quantity P 
and derive the order of the displacements as PL 2 j Eh from an asymptotical analysis 
of the energy functional. 

For a thin-walled slender beam whose dimensions satisfy Eq. (4.1) the rate of 
change of the displacements along the axial direction is much smaller than their 
rate of change along the circumferential direction. That is, for each displacement 
component 


dvi 

dx 


« 


dvi 

~d 7 


(4.15) 


Using Eq. (4.7) and assuming that d is smaller or of the same order as i£, the 
order of magnitude of the in-plane strains and curvatures is 


7 .. - 0 (|) 

2m ~ 0 (|) 
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Since 7 jj and p n are much smaller than 712 , 722 and pu, p 2 2 , respectively, their 
contribution to the elastic energy is neglected. 

The order of magnitude of the shell energy per unit area and the work done by 
external forces is 



Since PjUi << $, the contribution of external forces is neglected. Therefore the 
energy functional takes the form 

2 1 = J 0 ^{4/jC 1212 (7n) 2 ■+ 4/iC' 1222 7i2722 + hC 2 ~~ 2 [ ’ 122) 7 + 4/? 2 Cj 1212 7i2Pi2 

+ 2/j 2 C 1 1222 7 12 p 2 2 + 2ft 2 C 2212 722P 12 + h 7 C” 22 l22P22 
4 4 jcr’pneu 4 (4.16) 

Using Eq. (4.15), the strain-displacement relationships in Eq. (4.7) can be written 


as 


n dv l 

27,1 17 

dv 2 v 

7 ” = 17 + r 

1 dvy 

p " = 7117 

d 2 v d v 2 

Pz2 = l?~ls { R 


(4.17) 
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The integrand in Eq. (4.16) is a positive quadratic form, therefore the minimum 
of the functional is reached by functions v, t»i, and v 2 for which 712 = 722 = P12 = 


P 22 ~ 0. From Eq. (4.17) this corresponds to 


dv 2 v 

~d7+R- 

8a 2 ds\R) 

The function t> in Eqs. (4.19) and (4.20) should be single valued, i. e. 


( dv \ 1 f dv 

U) s 7 / Ts is = 0 


(4.18) 


(4.19) 


(4.20) 


(4.21) 


The bar in (4.21) and in the subsequent derivation denotes averaging along the closed 
contour T whose length is denoted by l in Eq. (4.21). 

Equation (4.18) implies that t»i is a function of x only. i.e. 


t'i = (M*) 


(4.22) 


Integrate Eq. (4.20) to get 


dv v 2 


(4.23) 


where <p(i) is an arbitrary function which is shown later to represent the cross- 
sectional twist. From Eq. (4.21) and (4.23), one obtains the relation between <p(x) 


and t> 2 . 


* (l) = (I) 


Substitute v from Eq. (4.19) into Eq. (4.23), to get the following second-order differ- 


ential equation for t>2 


d ( D dv 2 \ t'2 _ . 
6s ds j + R ^ 


(4.24) 
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To solve this equation, one has to recall the relations between the radius of curvature 

R and the components y(,s) and z(s) of the position vector associated with contour T 

<Pz 1 dy 

ds 2 R ds 

dPy 1 dz 

ds 2 R ds 


(4.25) 


It follows from Eqs. (4.25) and (4.4) that ^ and ^ are solutions of the homogeneous 
form of Eq. (4.24) and t> 2 = <p(x)r„ is its particular solution. The general solution is 
therefore given by 

t> 2 = U 2 (x + U 3 {x)^f- + y>(x)r„ (4.26) 

ds ds 

where U 2 and U 3 are arbitrary functions of x. Substitute from Eq. (4.26) into Eq. 
(4.19) to get 

v = U 2 (*)j^ ~ “ V(x) r t (4-27) 

Eqs. (4.22), (4.26) and (4.27) represent the curvilinear displacement field that mini- 
mizes the zeroth-order approximation of the shell energy. Using Eq. (4.8) the curvi- 
linear displacement field is written in Cartesian coordinates as 


vi = Ui(x) 


u 2 = U 2 (x) - z<p{x) 
«3 = U 3 {x) + y<p{x) 


(4.28) 


The variables Ui(x),U 2 {x) and U 3 (x) represent the average cross-sectional transla- 
tion while <p(x) the cross-sectional rotation normally referred to in beam theory as 
the torsional rotation. This displacement field corresponds to the zeroth-order ap- 
proximation and does not include bending behavior. For a centroidal coordinate 
system Ui(x), U 2 (x), U 3 (x) and can be expressed as 


U 2 (x) = ST 
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U 2 (t) = u 2 

U 3 (x) = til 

(v -t) 

V>(*) = ~~ 


(4.29) 


4.4.2 First-Order Approximation 

A first-order approximation can be constructed by rewriting the displacement field in 
Eqs. (4.22). (4.26) and (4.27) in the form 


Vi = Ui(x) + tl>i(i,x) 

v 2 = v 2 {x)~ + f/ 3 (*):r + v>(*)n, + ^'2(5,®) 

as as 

t» = - f ' 3 ( x )^7 ~ ¥ > ( a ’) r « + «’(•*, x) 


(4.30) 


where w } ,u> 2 and w can be regarded as correction functions to be determined based 
on their contributions to the energy functional. 

Substitute Eq. (4.30) into (4.7) to obtain the strains and curvatures in terms of 


the displacement corrections 


7n = Tn + 


o du> 2 dith 

2*)i2 = 27n 4- + 2-73.2 , 2712 = -7^- 

o . . dw 2 w 

722 = T22 + 722 » 722 = -qJ- + 


Pn = + gp- 

o d 2 u> 3 5 u> 2 . . 1 dtUi 

^>12 = P12 + ^-^7 - TZ~Zr + P** ’ 


5sdx 4i? 5x 


0 . # 2 u> 5 

Pn=Pn + hi , ** “-a?- Siva) 


where 7° a ^ and P° a $ are the strains and curvatures corresponding to the zeroth-order 
approximation. These are expressed as 


7n 



2 7„ = + <e\xy„ - o (|) 


722 = 0 






U' 2 ( X )^s + U ^ X ^7s + ^ X ) r " 



#22 — 0 


(4.32) 


The prime in Eq. (4.32) denotes differentiation with respect to x. Among the new 
terms introduced by the function u>, the leading ones are denoted by superscript * in 
Eq. (4.31). The order of ir, is (^). this is derived from Eqs. (4.31) and (4.32) where 
it is seen that the leading terms 2-712 and p J2 are of the same order of magnitude as 
2 7°i2 and P°i 2 , respectively, i.e. 


2712 — 


dwi 

1h 


. 1 d w i 

Pu = 7ri~d7 



(4.33) 


Therefore, 

n( Ad 

w ' ~ 0 (T 

An alternative approach is to assume the order of u>, as j and verify this assump- 
tion, as shown later, once iu, is determined. The order of magnitude of the remaining 
leading terms in Eq. (4.31) is as follows 


(4.34) 
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h '~' 0 {ji) (4 - 35) 

Tlie energy functional can be represented by $ ('>11,2712,722, Pu, p\2, 922)- By 
keeping the strains and curvature associated with the zeroth-order approximation 
and the leading terms contribution over the other terms (i.e., by dropping the terms 
and fe - in Eq. (4.31)) the energy function can be written as 

$0 11)21 12 + 2^12)0 + 722)^11)^12 4- Pl2)0 -f P22) 

The interaction terms associated with P°n and P°n, namely 

o 0 2 0 2 °** 

h P\\f >\2 > h P11P22 

o ^ °* 2 ° - 2 0 

hP 12712 9 hP 12I22 9 " P12P12 9 ^ P12P22 

are of order (757) or smaller. They are neglected in comparison with the following 
terms 

I'n7i2 ) I11722 ) I12712 ) I12722 (4.36) 

of order (jt). Similarly, the contribution of the work done by external forces, P,u,, is 
neglected since its order is (E/j jtC j)) in comparison with the order of the remaining 
terms in the energy functional ( Ehj y). Therefore in order to determine the functions 
Wi one has to minimize the functional 

j > $(l'ii)27i2 + 2712 ) 722) 0)Pl2)^22) d s 

If the rigid body motion is suppressed the solution is unique. The terms p 12 , P22 are 
essential to the uniqueness of the solution; however, their contribution to the energy, 
expressed by the interaction terms 

hpill 11 9 ^12^129 hpitflD ^P 22^12 


is of order ( 32 ( 5 )) or smaller, and is consequently dropped in comparison with the 
membrane contribution listed in (4.36). This aspect is discussed by Berdichevsky and 
Misiura [63], with regard to the accuracy of classical shell theory. Therefore, the shell 
energy can be represented by 


I — J j 1 + 2512 , 722 , 0 , 0 , 0 ) dsdi 


(4.37) 


It is worth noting that the bending contribution does not appear in Eq. (4.37). That 
is, to the first-order approximation the shell energy corresponds to a membrane stale. 
The first variation of the energy functional is 
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= ft/f 9* f (dvH \ 4 e*_ t (&», ^ «■ 
Jo J \d( 27 i 3 ) \ ds ) d^ 22 V ds 


+ R \ dsd ' T 


Recall that - a ^~^ = N i2 and ~~ = A r 2 2 , Eq. (4.38) takes the form 
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d{6wi) 


d(6u>2 ) 1 


5s 


Set the first variation of the energy to zero, to obtain the following 

dKu 


ds 

dN 2 2 
8s 
N22 
R 


= 0 

= 0 

= 0 


which result in 


and 


N 12 = constant 


(4.38) 


(4.39) 


N 22 — 0 


(4.40) 


This is similar to the classical solution of constant shear flow and vanishing hoop 
stress resultant. By setting N 22 to zero the energy density is expressed in terms of 
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7n and 71 2 only 

24>i = min 2$ = min hC* $yf ') a 0 0 ht = ^(^)(7n) 2 4 2 jB(s)7u7 12 4 C'Mtm)* (4.41) 

"Y22 T22 

The variables A{s), B(s) and C(s) represent the axial, coupling and shear stiffnesses, 
respectively. They are defined in terms of the D 0 ^ 6 as follows 

*»>=< D '"' > - ! <J^ >- ~ O(Eh) 

B(s) = 2 [< D" n ~ < D ' l ”£L D r’ > | - 0 ( EI>) (4.42) 

C(s) = 4 [< D m > > -O(Eh) 

where the 2D Young’s modulus D are expressed in terms of the Hookean tensor . 
E aB ’ yf ' in Eq. (4.10). The pointed brackets denote integration over the thickness as 
defined in Eq. (4.12). 

For convenience, D a ^ f is given in matrix form as 

\D] = [C? 1 ] - 2 [<?’] [<?’]"’ [C 4 ] +[<?"]’ [<?’] [<?-] (4.43) , 


W’here 


■D uu D u22 D uu 
\D] = D U22 D 2222 D 1222 

271112 pl222 £>1212 

$11 $12 $16 
[$ ] = $12 $22 $26 

.$16 $26 $66. 
$13 $15 $14 

[$ ] = Q 23 Q 2h Q 24 

.$36 $56 $46 J 


$33 $35 $36 

[$ 3 ] = $35 $55 $45 

. $36 $45 $44 J 


(4.44) 
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Ql3 Q 23 Q 36 
[^ ] = (? 15 Q 2i Q 66 

.Ql4 Qi4 Q 46 ■ 

The indices adopted in £q. (4.44) follow the convention of Ref. [50] . The bars over 
the reduced stiffness coefficients Q,j of Classical Laminate Theory, Refs. [19] and 
[50], indicate that these quantities are to be obtained through appropriate coordinate 
transformations. 

Equation (4.41) indicates that, to the first-order, the energy density function is 
independent of functions u> 2 and w. That is, the in-plane warping contribution to 
the shell energy is negligible. The function u>i however, can be determined from Eq. 
(4.39) and (4.41) and by enforcing the condition on to be single valued as follows 


1 

N 12 = -577 : = 7 (£(sbn + C , (s)'y 12 ) = constant 

0(2712) 2 


(4.45) 


Substitute from Eqs. (4.31) and (4.32) into (4.45) to get 




H 


UM 


dy 

ds 


dz 


dw 2 d\v\ 


+ + *•(*)>■.<*) 4 __+ ds 


constant (4.46) 


Following the relations in Eq. (4.15), the term ^ is neglected in comparison with 
7^. Moreover, the term in Eq. (4.46) may be neglected in comparison with 

This is possible, if [f?j is less or of the same order of magnitude as C. For 
the case when \B\» C additional investigation is needed. Since the elastic energy 
is positive definite, B 2 < AC, and B could be greater than C only if A >> C. In 
practical laminated composite designs \B\ < C as the shear stiffness is greater than 
the extension-shear coupling. Therefore, Eq. (4.46) becomes 


4 E/a(x)^ 4 ^'(*) r „( 5 ) 4 = constant 


(4.47) 
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Equation (4.47) is a first order ordinary differentia] equation in W\. The value 
of the constant in the right hand side of (4.47) can be found from the single value 
condition of the function uq : 



dw-i 

ds 


ds = 0 


(4.48) 


The solution of Eq. (4.47) is determined within an arbitrary function of x. This func- 
tion can be specified from various conditions. Each one yields a specific interpretation 
of the variable For example if tUf = 0 the variable Ui = t>7 according to Eq. (4.30). 
The choice of these conditions does not affect the final form of the ID beam theory 
and therefore will not be specified in this formulation. The result is the following 
simple analytical solution of Eq. (4.47) 


v>i = -yUfc) ~ - V-si*) + <?(s)*?'(x) + 9i( s Wi( x ) (4.49) 


where 

G(s) = j f [^y-c(r) - r »>( r ) dr ** 0 ( d 2 ) 

b 

r) c(r ) dr ~ O (d) 

c 

= '<’>=ck = ^ (4 ' 50) 

The area enclosed by contour T is denoted by A t in Eq. (4.50). It is seen from expres- 
sion (4.49) that Wj is of order and relation (4.34) is justified. The displacement 
field corresponding to the first correction is obtained by substituting Eq. (4.49) into 
Eq. (4.30) and dropping w 2 and w since their contribution to the shell energy is negli- 
gible compared to tnj . The result referred to as the first-order approximation is given 
by 

v, = V,(x) - y(s)Vfc) - z(t)UHx) + G(«)*>'(*) + «,(.)p;(») 

vi = 



(4.51) 
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v = U 2 {x)^f- - U s (r )^j- - <p(z)r, 
as as 


4.4.3 Second-Order Approximation 

Following a similar procedure to the one described in section 4.4.2, a second-order 
approximation can be constructed by rewriting the displacement field in Eq. (4.51) 
in the form 

Vi = Uj(t) -yU'(r) - zU'(x)-h G(s)tf'(z) -I- g 1 (s)U' 1 (x) + u'i(s.x) 


dy dz 

t» 2 = U 2 (x)— + ^(x)— + P( JP ) r « + tu 2 (« t ar) 
dz dy 

v = f 7 2(T)^ ~ U 3 ( T )fo ~ p( T ) r > + 


(4.52) 


where , u> 2 and w can be regarded as correction functions to be determined based 
on their contributions to the energy functional. 

Substitute Eq. (4.52) into (4.7) to obtain the strains and curvatures in terms of 
the displacement corrections 


7n = 7n 4 


du>i 

dx 


dw? - ■ dti^i 

2712 = 25l2 4 -TT~ 4 27 22 ? 27 2 2 = 


dx 


ds 


722 = 722 4 7 


22 1 


1 22 — 


dw? U> 
~ds + R 


Pn = pn 4 


d 7 w 

d ? 


(4.53) 


P12 — P12 4 


d 7 w 3 dw 


dsdx 4 R dx 


+ P 12 


P 12 = 


1 dw j 


4 R ds 


P22 — P22 4 p 22 , p 22 - 


d 7 w d (w 2 

Ts 7 


ds \ R J 


where and p a p are the strains and curvatures corresponding to the first-order 
approximation. These are expressed as 

JPL 


0-2. 
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2 W = ^cy'<*)+^t',V)~0(!) 


■722 = 0 

", Jy 


Pn = W(x)f s - UZ(x)-jL - ¥>V)r, ~ o (-) 


(4.54) 


P22 — 0 

The terms written over the overbraced expressions in Eq. (4.54) denote their order. 
Among the new terms introduced by the function th, the leading ones are denoted by 
superscript" in Eq. (4.53). The order of th, is assumed to be 


th, ~ O 


ffl 


(4.55) 


Such an assumption will be justified later. Therefore, the order of magnitude of the 
leading terms, Eq. (4.53), is as follows 


: : n (Ad 

7 12 ~ 722 ^ 0 I — 


Pl2 ~ ^22 ~ 0 


(4.56) 


The energy functional can be represented by $ {‘jn->‘Zli 2 ,l 22 iPu,Pi 2 ,P 22 )- Bv 
keeping the strains and curvature associated with the first-order approximation and 
the leading terms contribution over the other terms (i.e., by dropping the terms 
|p, and i n Bq. (4.53)) the energy function can be written as 


$(7lli27i2 + 27 i 2 ,0 -f 722>frl>Pl2 + Pi2;0 + P 22 ) 


In the following, the order of magnitude of the energy due to bending, i.e. due to pu, 
P 12 , P 121 p 2 2 , is investigated. 
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The interaction terms associated with pn, namely 


"711712 i 7ll7"22 5 7l27l2 i 712722 i 


(4.58) 


^Pn7i2 , hpij7'j2 5 ^ P\\Pi7'> h P 11 P 22 

are of order (^r^) or smaller. They are neglected in comparison with the following 
membrane contribution to the energy 

0 (^) associated with V[ and p' 

0 associated with U 2 and U£ 

The interaction terms due to the bending curvature p J2 are 

hpi 27 i 2 » hp 12 i 22 ~ O associated with U[ and p' 

h 2 pi 2 Pj 2 ? h 2 p 12 p 22 ~ O (^)] assoc ’ a<ec * ^ (4.60) 

These terms are of higher order in comparison with the membrane contribution asso- 
ciated with U[ and ip' in Eq. (4.58), and may be neglected. The remaining interaction 
terms associated with p 12 and p 22 . namely 

O ( 7 T 1 ) associated with U[ and p' 


(4.59) 


hhiPii j h 7 jip 22 , h'fi 2P12 1 ^712^22 1 


rr" 

3 


(4.61) 


[ ~ 0 associated with U 2 and U 

are of higher order when compared to the corresponding membrane ones, listed in 
(4.58). Therefore in order to determine the functions u.*, one has to minimize the shell 
energy expressed by 

rL 


I = J o j $( 711*2712 + 2 -) 12 , 7 22 , 0 , 0 , 0 )dsdx 


(4.62) 


The contribution of the new corrections in the work done by external forces is neg- 
ligible compared to the first-order approximation. Consequently its contribution is 
neglected in Eq. (4.62). It is worth- noting that the bending contribution does not 
appear in Eq. (4.62). That is, to the second order approximation the shell energy 
corresponds to a membrane state. 
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The first variation of the energy functional is 

«/= f l l{ g* f'li 

Jo J 15(2,,:) \ds ) dm 1, 3s + Rjf 

Recall that y = Ny 2 and = N 22 , Eq. (4.63) takes the form 
cr t L I \ * r d(6w!) , XT (d(6w 2 ) . 1 ,. r M , 


(4.63) 






Set the first variation of the energy to zero, to obtain Eqs. (4.39) and (4.40). By 
setting A t 2 2 to zero, the energy density is expressed in terms of 711 and -) 12 only as 
given by Eq. (4.41). The function can be determined from Eq. (4.39) and (4.41) 
and by enforcing the condition on ibi to be single valued as previously outlined in 
section 4.4.2. Substitute from Eqs. (4.53) and (4.54) into (4.45) to get 

(a£V| 

(EH) (ff) Iff) 

\ 'b' VUx)- - z{,s)iq 4 G(s)?"(x) + S,( S )P;'(I) + 


<*> (#) (£)’ 


1 2 >l ( , dpi , c?u» 2 5u‘i 

+ 4 C 


= constant 


(4.64) 


Comparing the order of magnitude of each kinematical variable, Eq. (4.64) reduces 
to 


+ \c -£rcp'(z) + = constant (4.65) 

4 ic as os 

Using the single value condition of function tbi, the following simple analytical solution 
of Eq. (4.65) is obtained 


£>i=92UZ(x) + g 3 (s)UZ(x) 


(4.66) 
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where 

p 2 (js) = -jT i(r)y(r) - y c ( t ) ~ ^ (^ 2 ) 

S3( s ) = ~ J q h(T)r(r) - y c ( t ) ^ 0 (rf 2 ) (4.67) 

It is seen from expression (4.50) and (4.67) that G{s). gi(s), 32(5) and s 3 (s) are 

single- valued functions, with 

G( 0 ) = G(l) = pi(O) = 31 (I) = 32(0) = g 2 (l) = 33(0) = g 3 (l) = 0 

Using Eqs. (4.66) and (4.67), tl>i is found to be of order and the assumption in 

Eq. (4.55) is justified. 

4.4.4 Convergence of Displacement Field 

The displacement field corresponding to the second correction is obtained by substi- 
tuting Eq. (4.66) into Eq. (4.52) and dropping u> 2 and u; since their contribution to 
the shell energy is negligible compared to tt»i. The result is 

t’l = EM?) — V(s)Vl(x) - z(s)U'(x) -f g(s)yV) 

+ 9i(*)Fl( x ) + 32(^)f r 2 , ( ;r ) +9a(»Ws( r ) 

t'2 = + r/ 3 (x)~ + ¥>( r K (4.68) 

as as 

v = U 2 {x)-^-U 3 (x)^-ip(x)r t 

A third cycle is carried out, however no additional terms of the same order in the 
energy functional result as shown in the Appendix, and the final displacement field 
converges to the expression given in Eq. (4.68). 

The underlined terms in Eq. (4.68) correspond to the classical theory of extension, 
bending and torsion of beams. The additional terms gi(s)U[, giis)^' and g 3 {s)U 3 
in Eq. (4.68) represent warping due to axial strain and bending . These new terms 
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emerge naturally in addition to the classical torsional related warping G(s)<p'. They 
are strongly influenced by the material anisotropy, and vanish for materials that are 
either orthotropic or whose properties are antisymmetric relative to the shell middle 
surface. For these layups the coupling parameter h(s) defined in Eqs. (4.50) and 
(4.42) vanishes. The significance of the axial and bending-related warping terms and 
their effect on the accuracy, is shown in the applications of Chapter V. Moreover, the 
expression for torsional related warping G(s) differs from the work of Refs. [30] and 
[42]- [46] . A comparison of these expressions is presented in section 4.6. 


4.4.5 Strain Field 


We now substitute the displacement field of Eq. (4.68) into the in-plane strain com- 
ponents of Eq. (4.7), while using Eq. (4.50), to obtain 


JilL ML (fj!) 

*>11 = Ui(*)- yUt(x)- zUH(x)-\- G(s)<p"(x) 

m {¥■) Jg) _ 

+ 9i(*W"(x) + 92V?(z) + gtUS’ '(*) 


(£) 


(f) 


(**) 


i&) 


2li2 = 


// V , ^9 1 rr/ / v , ^2 rf// / \ , jrlt f V 

—cr(z)+—v^) + —l,(x) + —lAT) 


dgi 


dg2 


dg 3 


722 — 0 


(4.69) 


The terms giU", giU™ , and g 3 Uz can be neglected in comparison with U[, yU'i, and 
zU$, respectively. Therefore, the in-plane strain components become 


7» = U[(x) - y(s)U^x) - z(s)U!(x) + Gtp" 


o 1A t , dg\ , dgz TT „ dg$ „ 
2-m = -JfV + ^ + -^V, 4 — U 3 


(4.70) 


722 — 0 
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Using Eq. (4.70), the shell energy density, Eq. (4.41), can be written as 
(EM 

2#! = (U[f + (yU''f 4 (zU‘‘f -f yU'JJZ - 2zV[U« 

■ m j . 

+ - 2yGV& - 2:GI/,y 


(*) 




: -37 (ir »i +G ir c ^ 9 


i ( i ^?TT'TT"i ^93 T jt TT n 2A f u , 

+ -jj-cT'.io +.-57^.K, + 

- T T "V' — d'9& TT"TT" ft" r! . ^91 TT tt i t/ 

s, 77 l ^ ' v j7 Ui v * ~ TT c - £ 3 10 - ! 3 ' ■ 

(«*) (w m 1 

■ s /■ ■ ■ | A " - 'S /"■ s 

fjtifj/t I Jfrr 1 . r* **'92 MtjH , /~i ^9z //tt// 

-■-— +G— y P.+g^ ^ 


^ <M*(H + (**)’♦ (M 


(#) 


(#) 


4A e rr/ , , 4A e dp 2 rr „ , 4A e dp 3 rr „ , 

4 4 TT^’ 1 ’’ 4 lf c ^ u ^ 

4 2 1717 u M 4 *££«Wr + 2%%y"£/; 


ds ds 


ds ds 


(4.71) 


where the underlined terms are associated with the GV> ,/ contribution in Eq. (4.70). 
These terms are of higher order and may be neglected in comparison with the remain- 
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mg overbraced leading terms, as shown in Eq. (4.71). Therefore, one may drop Gf" 
from Eq. (4.70), and the final expressions for the in-plane strain components, using 
Eqs. (4.50) and (4.67), become 

lu ^ U[(x) - y(s)U^x) - x(s)U^x) 


2*)12 = 


2-4-e , \ / , 

— c(s)^ + 


6(s) - = c(s) 
c 


v[ 


b{s)y(s) ~ ~ c i s ) 

c 

6(s)r:(s) - ^c(s) 
c 


V" 
1 2 


722 = 0 


(4.72) 


It is worth noting that the vanishing of hoop stress resultant in Eq. (4.40) and hoop 
strain in Eq. (4.72) should be interpreted as negligible contribution relative to other 
parameters. The longitudinal strain 7 u is a linear function of y and z. This result 
was adopted as an assumption in the work of Ref. [29]. 

In the present formulation, parameters A , B and C where assumed to be of 
the same order. However, the results are valid for configurations which satisfy the 
following inequalities 


A 

C 



« 1 



<< 1 


4.4.6 Constitutive Relationships 


Dropping the underlined terms in Eq. (4.71) and integrating over the shell middle 
surface to get the energy of ID beam theory 


I 





-/ 


PiUidxd-s 


(4.73) 
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where 

= “ [C'ii(^j) J -f C 22 (<p'f + C 3 3(Us)‘ + C'44([/j ) 2 ] 

+C’ 12 ^V + c l3 uw + C 2 A U[IV' 

+c 3 ztp'UZ + CWt? + CWOT (4.74) 

Explicit expressions for the stiffness coefficients C,j (i, j = 1, 4) are given in Eq. 
(4.78). 

The constitutive relationships can be written in terms of stress resultants and kine- 
matic variables by differentiating Eq. (4.74) with respect to the associated kinematic 
variable or by relating the traction T , torsional moment M x . and bending moments 
M y and M z to the shear flow and axial stress as follows 


T =W{ = f!‘ ,ndids =l' 

' N u ds 


M x = d , = / J (*)d£ds = j 

£ Ni 2 r„(s)ds 


M ' = Ws = = - 

j> Nnz(s)ds 

(4-75) 

d$ 2 r r 

M z = ^777 = ~ f J vuydids = - 

j Nuy(s)ds 



The shear flow N i2 is derived from the energy density in Eq. (4.45) and the axial 
stress resultant N u is given by 

Ah = = A(s)lu + B(s)~n 2 (4.76) 

07 u 

and the associated axial and shear stresses are uniform through the wall thickness. 
Substitute Eq. (4.72) into Eqs. (4.45) and (4.76) and use Eq. (4.75) to get 


' rp ’ 


Cu 

Cl 2 

C13 

C14 



M. 

, _ 

C i2 

c 22 

C 2 3 

C 2 4 

< 


My 


C13 

C 23 

C33 

C34 



\M Z \ 


Cl 4 

C 2 4 

C34 

C44 _ 


l^'J 


(4.77) 
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where expressions for the stiffness coefficients Cij (i, j — 1, 4) in terms of the cross 
section geometrj' and materials properties are as follows 

due to g] V[ 


C-n = / M - ~)ds 4 


B\, , | f(BIC)dsf 

~C ] 4 7(17 C)Z 


" aiic)is A ‘ 


due to g\V[ and gsVy 


— — j>(A — ^r)zds ~ 


f(B/C)ds f(B/C)zds 
/(I /C)ds 


due to gi V; and gjV" 


Cj 4 = - j>( A - Yj-)yds - 


§(B/C)ds $[B/C)yds 
f{l/C)ds 


C 22 f(llC)d* A * 


= 


f(B/C)zds 
fd/C)ds ' 


/ (B/C)yds 4 

Cu ~ -~IWW A ' 


C 33 = £{A 4 


due to g}Uy 

~ TfWC):isf 

#(1/0* 


C34 — j>(A — — )yzds ■+ 


due to g^Uj' and gjU^ 

B r , . , JIbJc^Is TTb/ct^ 


f(l/C)is 


due to gjVj 


c « = f ( A ~ ¥r)y 7 ds + 


l. §(B/C)yds ] 2 


(4.78) 


" J v C’ S i f(l/C)ds 

The out-of-plane warping contribution to the stiffnesses due to the axial strain (i.e., 
due to g\U[), bending about y axis (i.e., due to gzU^), and bending about z axis (i.e., 
due to < 72 ^ 2 ) is shown by the overbraces in Eq. (4.78). 
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The coefficients C tJ (i, j = 1,4) can be expressed in terms of the in-plane axial 


stiffness coefficients yl tJ of Classical Lamination Theory (CLT) if one neglects the 
through-the-thickness contribution to the stiffnesses in Eq. (4.78). The result is 


c u-f(h n -—)is + H i jKnVtt _ 


0 ^ ( Ai 2 /A 22 )ds 

C ' 12 2>le / (l/A '22 )ds 

„ A'* lv , §{KJK 22 )ds §{KJK„)zds 

t -13 = - f (All - -zr~) zds 777 i r " j 
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where, the stiffnesses A\j are 


An — -^n ~ 


(^ 12 ) 


A 12 — Ai6 — 


j 4 i 2 j 4 2 


(4.79) 


K 22 = A 6 6 — 


(■'426 ) 


4.4.7 Equilibrium Equations 


The equilibrium equations are derived using the principle of virtual work. The vari- 
ation of the internal strain energy is 

61 1 = J o f'(Kzxhxx + 2N x ,hx.) dsdT 

Using the strain displacement relations, one-dimensional stretching, twisting, and 
bending generalized internal forces are defined as 

T = f N xx ds 

M x = 2A t N x , 

M y = — £ N xx zds 
M : — - j> N xx yds 

Consider a beam subjected to external forces and moment resultants T, ~M X , M y . and 
M. at both ends. Moreover, surface tractions P x , P y , and P z are applied along the 
i, y, and z directions, respectively. The variation of the virtual work of the external 
forces can be written as 

m = -My8Uti-M z 6Uti 

4- £ [(^ P x dsj 6U 1 - (y P x yds ) 6(7' - ^ P x zds^ 6U' 3 + P v ds ) 6U 2 

— [j> Pyzds^j 6ip — P z ds j 6U 3 — P z yds j 6pj dx 

Using the principle of virtual work 

6U = 6W e . 


one obtains a system of linear equilibrium equations as follows 
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M' r + f (P z y - P y z)ds = 0 

A/; -f {j P r zds)' + f P z ds = 0 (4.80) 

m; + ^ p v ds = o 

One of the member of each of the following four pairs must be prescribed at the 
beam ends : 

T or Vi , M x or ip, M y or U^, and M z or V' 2 (4.81 ) 

4.5 Summary of governing equations 

The development presented in this work encompasses five equations. The first, is 
the displacement field given in Eq. (4.68). Its functional form was determined based 
on an asymptotical expansion of shell energy. The associated strain field is given in 
Eq. (4.72) and the stress resultants in Eq. (4.45), (4.75) and (4.76). The fourth, are 
the constitutive relationships in Eq. (4.77) with the stiffness coefficients expressed as 
integral of material properties and cross sectional geometry in Eq. (4.78). Finally the 
equilibrium equations and boundary conditions are given in Eqs. (4.80) and (4.81), 
respectively. 

In the present development the determination of the displacement field is essential 
in obtaining accurate expressions for the beam stiffnesses. A comparison of the derived 
displacement field with results obtained by previous investigators is presented in the 
following section. 

4.6 Analytical comparison with previous results 

In anisotropic materials the importance of physical effects such as transverse shear 
strains is influenced by the relative magnitude of elastic moduli. For example in 
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laminated composites the extensional modulus along the fibers direction is usually 
large relative to the shear moduli and consequently transverse shear effects can be 
significant. Several theories have addressed this issue by including transverse shear 
in the assumed displacement field [30], and [42]- [46] . The displacement function 
Eq. (4.68) derived from the asymptotical analysis does not include transverse shear 
strain terms explicitly. This is a consequence of the vanishing of the through-the- 
thickness stress component <r* 3 in Eq. (4.5) or (4.9) where the transverse shear strains 
are expressed in terms of other strain components. Their effect however is implicitly 
included in the stretching-related warping term gi{s) and the bending- related warping 
terms p 2 (.s) and 53(5) as illustrated by the applications of Chapter V. 

Rehfield’s theory [30] recognizes the significance of transverse shear strain in thin- 
walled composite beams. Its displacement field is given by 

«i = £M*) -y(s)[f' r 2 (*) - 27* v (t)] - z(s)[Us{x) - 27 ,-(i)] + g(s,x) 

«2 = l r 2 {r) - z{s)tp{x) (4-82) 

U3 = U 3 {x) + y(s)ip(x) 

where 7^ and 7, ; are the transverse shear strains. The warping function g(s, x) is 
given as 


g{s,x) = G{s)<f'(x) 

(4.83) 

G(s) = 2A e y- /*r„(r)(fr 
l Jo 

(4.84) 


A comparison of the displacement fields in Eq. (4.68) and (4.82) shows that the 
warping function in Rehfield’s formulation includes the torsional-related contribution 
and does not include explicit terms that express the bending-related warping. The 
torsional-related warping function G(s) in Eq. (4.50) is different from the function in 
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Eq. (4.84). The two expressions coincide when c = constant that is, when the wall 
stiffness and thickness are uniform along the cross section circumference. 

The torsional related warping function in Eq. (4.84) was modified by Atilgan [44], 
and Rehfield and Atilgan [43] as 


f » TO A 

G{s) =L fe c, - r " (T) . 


dr 


(4.85) 
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(4.87) 


The A{j in Eq. (4.87) are the in-plane axial stiffnesses of CLT, Refs. [19] and [50], 
they are related to the modulus tensor by 


A u =< jE 1111 > , .4 12 =< E u22 > , A 22 =< E 2222 > 

Aie =< E ul2 > , A 26 =<E 1222 > , A 6 s=<£ m2 > 

A comparison of the modified torsional warping function in Eq. (4.85) and G{s) in 
Eq. (4.50) shows that they coincide for laminates with no extension-shear coupling 
( < D im > = < D 1222 >= 0, in Eq. (4.10) ). For the case where the through-the- 
thickness contribution is neglected in Eq. (4.10), this reduces to Ai 6 = A 2 e = 0. 

The warping function obtained in Refs. [42] and [46] for composite box beams is 
identical to the expression of Refs. [43] and [44] in Eqs. (4.83) and (4.85). 

An assessment of all the previous warping expressions can be made by checking 
whether they reduce to the exact expression for isotropic materials (see, for example, 
Ref. [59]) 

< T >h 


(4.88) 


where // is the shear modulus. 

For isotropic materials the in-plane coupling b is zero and consequently g 2 and 
gz in Eqs. (4.50) and (4.67) vanish. That is the warping is torsi on -related and reduces 
to G(s)ip'. Moreover, the shear parameter c is equal to and the expressions for 

(?(.s) and G(s) in Eqs. (4.50) and (4.88) coincide. 

Rehfield's warping function in Eq. (4.84) coincides with Eq. (4.88) when the ma- 
terial is isotropic and the wall thickness is constant. Also the works of Refs. [43], [44] 
and [46] reduce to Eq. (4.88) for isotropic materials. 

4-7 Closing Remarks 

The major advantage of the approach adopted in this work is the fact that the dis- 
placement function emerges as a result of the asymptotical analysis of the shell energy. 
The influence of the material's anisotropy is accounted for in a consistent manner and 
the deformation modes are determined on the basis of their contribution to the asso- 
ciated energ)\ Two new contributions to the warping emerge due to stretching and 
bending. They are of the same order of the classical torsional-related warping. Their 
significance is illustrated in the applications provided in the next chapter. 


97 


CHAPTER V 

APPLICATIONS OF ANISOTROPIC THIN-WALLED 

BEAM THEORY 


An evaluation of the variation ally consistent theory developed in chapter IV is 
provided. The theory is applied to beams with arbitrary closed cross-sections made 
of laminated composite materials with variable thickness and stiffness subjected to 
axial load, torsion and bending. A comparison of flexibility coefficients and deforma- 
tion with finite element predictions, closed form solutions and experimental data is 
performed to validate predictions and isolate the influence of different contributions to 
the section warping. In addition to the torsional related warping, two new contribu- 
tions namely, axial strain and bending related out-of-plane warping were identified in 
the developed theory. Extension and bending related out-of-plane warping are shown 
to have a significant effect on the accuracy of predictions. Comparison of predictions 
provides also a check of the asymptotical analysis result regarding the contribution 
of shear deformation. Although the resulting displacement field does not include 
an explicit shear deformation term similar to Timoshenko’s theory, shear deforma- 
tion contribution is shown to be implicitly accounted for through the out-of-plane 
warping due to extension and bending. 

Two special layups: The circumferentially uniform stiffness (CUS) and circum- 
ferentially Asymmetric stiffness (CAS) have been considered in Refs. [41]- [46] and 
[51]. They are associated with different non-classical behaviors. These behaviors are 
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shown to be influenced by the out-of-plane warping due to extension and bending in 
the next section. 

5.1 Effect of Out-of-Plane Warping due to Extension and Bending 


5.1.1 CUS Configuration 


This configuration produces both extension-twist and bending-transverse shear cou- 
plings. The axial, coupling and in-plane stiffnesses A , B, and C given in Eq. (4.42) 
are constant throughout the cross section and hence the name circumferentially uni- 
form stiffness (CUS) adopted in Ref. [43], [44], [45] and [51]. Such a configuration 
is manufactured by wrapping the composite lay-up using a winding technique. For 
a box-beam, the ply lay-ups on opposite sides are of reversed orientation, and hence 
the name antisymmetric configuration adopted in Refs. [41], [42], and [46]. 

Since A , B, and C are constants, the stiffness matrix in Eq. (4.78), for a centroidal 
coordinate system, reduces to 


Cji C'i 2 0 0 

C12 C'j 2 0 0 

0 0 < 7 33 0 

0 0 0 C44 


(5.1) 


The nonzero stiffness coefficients are given by 


C u = Al 


Ci 2 = BA e 
C22 — jA 2 t 

C33 = A £ z 2 ds — — z 7 ds 
r B 2 r 

C44 = A j y 2 ds - — j> y 2 ds 


(5.2) 
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where l denotes the length of the closed contour T. For such a case the out-of-plane 
warping due to axial strain vanishes and gi does not affect the response. This is 
shown by considering A , B , and C as constants in Eq. (4.78). The influence of the 
out-of-plane warping due to bending in the x-z and x-y planes are expressed by the 
underlined terms in the expressions of C 33 and C44, respectively. These terms are 
significant in predicting the deflection of antisymmetric configurations. 


5.1.2 CAS Configuration 


This configuration produces both bending-twist and extension-transverse shear cou- 
plings. The stiffness A is constant throughout the cross section. For a box beam, the 
coupling stiffness, B , vanishes for the vertical members, while its values in the top 
and bottom members are of opposite signs 

B top — “ B bottom 


Bvertieal members = 0 (5.3) 

and hence the name circumferentially asymmetric stiffness (CAS) adopted in Ref. [43], 
[44], [45] and [51]. For a box-beam, the ply lay-ups on opposite sides are mirror images, 
and hence the name symmetric configuration adopted in Ref. [41], [42], and [46]. The 
stiffness C along the horizontal and vertical members are equal and expressed by 


C top — C bottom 


C vertical left — C vertical r 


ight 


(5.4) 


The stiffness matrix, for a centroidal system of axes, reduces to 

Cu 0 0 0 


[Cij] = 


0 C22 0 23 0 

0 C 23 O 33 0 

0 0 0 C44 


(5.5) 
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The nonzero stiffness coefficients are expressed by 


Cu = Al-2^d 


(?22 ~ 
C'23 = 


C t 


»!*'+-(&)] ' 

h A * 

*'+■(£)] ' 



a 


44 


- A f> 


2 ds 


Bfd 3 
6 C t 


(5.6) 


Subscripts t and v denote top and vertical members, respectively. The box width 
and height are represented by d and a, respectively. Equations (5.6) are derived by 
substituting Eqs. (5.3) and (5.4) into Eq. (4.78) and considering A to be constant. The 
underlined term in the expression of the axial stiffness Cu represents the extension 
contribution to the out-of-plane warping. The bending contributions to the out-of- 
plane warping are represented by the underlined terms in the expressions of C33 and 
C'44 ■ For the CAS configuration, bending about the t/-axis is coupled with torsion 
while extension and bending about the r-axis are decoupled. 

In order to assess the accuracy of the predictions and isolate the influence of 
stretching and bending-related warping, the present theory is applied to the box 
beam given in Ref. [51]. The cross sectional configuration is shown in Fig. 5.1 and 
the material properties in Table 5.1. 


5.2 Comparison of Flexibility Coefficients 

A comparison of the flexibility coefficients Sij with the predictions from two models 
is provided in Table 5.2. The flexibility coefficients 5,j are obtained by inverting 
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Table 5.1: Properties of T300/5208 Graphite/Epoxy 


En = 21.3 Msi 
E 22 = E 33 = 1.6 Msi 
(?i 2 = G12 — 0.9 Msi 
G23 = 0. / Msi 
V\2 — 1^13 = 0.28 
i / 2 3 = 0.5 


the 4x4 matrix in Eq. (4.77). NABSA (Nonhomogeneous Anisotropic Beam Section 
Analysis) is a finite element model based on an ext ension of the work presented in Ref. 
[32]. In this model all possible types of warping are accounted for. The TAIL model 
is based on Ref. [30], but neglecting the restrained torsional warping. The predictions 
of the NABSA and TAIL models are provided in Ref. [51]. The percentage differences 
appearing in Table 5.2 are relative to the NABSA predictions. The present theory is 
in good agreement with NABSA. Its predictions show a difference ranging from +0.7 
to +3.6 percent while those based on Ref. [30] range from +3.6 to —18.4 percent. 

Since the box beam has a CUS configuration, the out-of-plane warping due to 
bending has a significant effect on the prediction of the bending flexibilities and 
(^) as shown in Eq. (5.2). Neglecting <73 and <? 2 in the expressions of C33 and C 44 
leads to values of 0.11424 x 10 -4 lb - 1 in -2 and 0.38410 x 10 ~ 4 lb _ 1 in _J for S 33 and 
5 44 , respectively. Comparison of these values with the underlined results in Table 5.2 
shows a 65 percent increase in the bending flexibilities due to out-of-plane bending 
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Table 5.2: Comparison of Flcvcibiliiy Coefficients of NABSA, TAIL and Present 
(lb, in units) 


Flexibility 

NABSA 

PRESENT 

% Diff. 

TAIL 

% Diff. 

S n x 10 5 

0.143883 

0.14491 

+0.7 

0.14491 

+0.7 

S 2 2 X 10 4 

0.312145 

0.32364 

+3.6 

0.32364 

+3.6 

S n x 10 5 

-0.417841 

-0.43010 

+2.9 

-0.43010 

+2.9 

S 33 x 10 4 

0.183684 

0.1886 

+2.6 

0.17294 

-5.8 

S 44 x 10 5 

0.614311 

0.63429 

+3.2 

0.50157 

-18.4 


related warping. 

5.3 Comparison of Deformation 

The present theory is applied to the prediction of the tip deformation in a cantilevered 
beam made of Graphite/Epoxy and subjected to different types of loading. The beam 
has a CUS square cross section with [+12]* lay-up. The geometry and mechanical 
properties are given in Table 5.3. Comparison of results with the MSC/NASTRAN 
finite element analysis of Ref. [38] is provided in Table 5.4. The applied axial and 
transverse forces are equal to 100 lb, while the applied torsional moment is 100 lb-in. 

The MSC/NASTRAN analysis is based on a 2D plate model accounting for both 
shear deformation and warping. The predictions of the present theory range from 
+1.7 to —0.7 percent difference relative to the finite element results. 

The deflection due to transverse load neglecting out-of-plane bending related warp- 
ing is equal to 1.341 inch compared to 1.853 inch (38% difference) in Table 5.4. For 
a CUS configuration, the extension-torsional response is decoupled from bending as 
shown in Eq. (5.2). Since C is constant and does not affect the stiffness coefficients, 



103 


Table 5.3: Geometry and Mechanical Properties of Thin-Walled Beam with [•+ 12] 4 
CUS square cross-section 


Length = 24.0 in. E u = E 22 = £ 33 = 11.65 Msi 

Width = depth = 1.17 in. G 12 = G 13 = 0.82, G 2 3 = 0.7 Msi 
Ply thickness = 0.0075 in. v 22 = u is — 0.05, ^23 = 0.3 


as outlined in section 5.1.1, the flexibility coefficients controlling extension and twist 
response, 5 ji, 5n and 522 coincide with those of Refs. [43] and [44]. As a conse- 
quence, the axial displacement and twist angle predictions coincide. However, the 
lateral deflection under transverse load differs. The tip lateral deflection predicted 
using the theory of Ref. [30], which includes shear deformation, and Refs. [43] and 
[44], which include a shear deformation correction to Ref. [30], is 1.724 inch resulting 
in —7.6 percentage difference compared to the NASTRAN result. This is due to the 
effect of bending-related out-of-plane warping on the bending flexibilities ^ and 
(O33 = C44 for this case), as shown by the underlined terms in Eq. (5.2). 

Figures 5.2 and 5.3 show the bending slope variation along the beam span for 
antisymmetric and symmetric cantilevers under a 1 lb transverse tip load, respec- 
tively. The beam geometry and its material properties are given in Table 5.5. The 
experimental results are reported in Refs. [41], [42], and [46]. The influence of the 
out-of-plane warping due to bending is isolated in these figures. The bending related 
out-of-plane warping, g^U 2 and gzU 2 terms in Eq. (4.68), results in a 91 and 20 % 
increase in the bending slope for the antisymmetric and symmetric configurations, 
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Table 5.4: MSC7NASTRAN and Present Solutions for a CUS Cantilevered Beam 
with [+12] 4 Layups Subjected to Various Tip Load Cases 


Tip Load 

Tip Deformation 

% Diff. 

Axial Force 
Axial Force 
Torsional Moment 
Transverse Force 

NASTRAN Present 
Axial Disp. : 0.002189 in. 0.002202 in. 
Twist : 0.3178 deg. 0.32325 deg. 

Twist : 2.959 deg. 2.998 deg. 

Deflection : 1.866 in. 1.853 in. 

+0.6 % 
+ 1.7 % 
+ 1.32 % 
-0.7 % 


respectively. The analytical predictions reported in Refs. [41], [42], and [46] together 
with results obtained on the basis of the analyses in Ref. [30], [43], [44] and the present 
theory are combined in Figs. 5.4 and 5.5. Results show that the present theory is 
in good agreement with the test data and the closest when compared to the other 
analytical approaches which include shear deformation, Refs. [30], [42], and [46], and* 
shear deformation corrections, Refs. [43] and [44]. 

The bending slope in Figs. 5. 2-5. 5 is defined in terms of the cross section rotation 
for theories including shear deformation. For the geometry and material properties 
considered, this effect is negligible as shown in Figs. 5.4 and 5.5 where the spanwise 
slope at the fixed end from theories with shear deformation, is indistinguishable from 
zero. The nonzero value shown by the test data may be due to the experimental set 
up used to achieve clamped end conditions. 

The spanwise twist distribution of symmetric cantilevered beam with j30] 6 and 
[45]6 lay-ups is plotted in Figs. 5.6 and 5.7, respectively. The beams are subjected to 
a transverse tip load of 1 lb. Their dimensions and material properties are given in 
Table 5.5. Results show that the present theory and those of Refs. [43] and [44] are 
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Table 5.5: Cantilever Geometry and Properties 


Width = 0.953 in. E n = 20.59 Msi, E 22 = £33 = 1.42 Msi 

Depth = 0.53 in. G J2 = G 13 = 0.87 Msi, G ' 23 = 0.7 Msi 

Ply thickness = 0.005 in. i / 22 = ^13 = 0.42, t / 23 = 0.5 


the closest to the test data. A similar behavior is found for the bending slope and the 
twist angle at the mid-span of the symmetric cantilevered beams appearing in Figs. 
5.8 and 5.9. The beams are subjected to a tip torque of 1 lb-in. 

5.4 Shear Deformation Contribution 

The significance of the out-of-plane warping due to bending is illustrated in Fig. 5.2. 
A similar behavior is obtained in Ref. [65] when the shear deformation contribution 
is neglected. This indicates that the out-of-plane warping due to bending includes 
implicitly the shear deformation contribution. In order to assess this similarity, the 
present theory and the numerical work of Ref. [65] are applied to the prediction 
of the deflection in a cantilevered beam made of grap hite/epoxy and subjected to 
a transverse tip load of 1 lb. The beam has a CUS cross-section with [+15]6 lay- 
up. The geometry and mechanical property, provided in Ref. [65], are given for 
convenience in Table 5.6. Figure 5.10 shows a similar behavior suggesting that in the 
present theory, shear deformation is implicitly accounted through ben ding- related 
warping. The prediction of Ref. [65] are referred to Classical when shear deformation 
is neglected. Further evidence could be provided by estimating the equivalent shear 
deformation strain. This can be expressed by the slope of the plane that approximat es 
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Table 5.6: Cantilever Geometry and Properties 


Width = 0.923 in. = 20.6 Msi, £22 = £33 = 1-42 Msi 

Depth = 0.50 in. G 12 = G 13 = 0.87 Msi, G 23 = 0.696 Msi 

Ply thickness = 0.005 in. 1/12 = t/u = 0.3, 1/23 = 0.34 


the cross-section warping and is given (66] by 

/ y Vi dA 


^■Hxy — 


L. 


(5.7) 


where A and I zz represent the cross-sectional area and moment of inertia about the 
r-axis, respectively. 

For a CUS box cross-section subjected to a vertical tip transverse load p., the 
shear strain distribution across the cantilever length is obtained by substituting the 
axial displacement t>i from Eq. (4.68) into Eq. (5.7). The result is the following 
analytical expression 

(L — * 1 ) had f a 2 . d 2 \ 

2hv = (y + + yl ^33 Pz (5.8) 

where 


S33 = Bending flexibility 
L = Length of cantilever 

X\ — Cross-section position measured from the fixed end 
h = Laminate thickness 
a = Box height 
d = Box width 
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A comparison of the shear strain 7 xy over the length of the cantilever with the 
prediction of Ref. [65] is shown in Fig. 5.11. The shear strain at the fixed end is 
4.5924 x 10 ~ 4 based on Eq. (5.8) which is within 2 percent of 4.6857 x 10 -4 calculated 
on the basis of Ref. [65]. 

5.5 Conclusion 

The anisotropic thin-walled closed section has been validated by comparison of re- 
sponse predictions with finite element solutions, other closed form analyses and test 
data. The influence of the two new nonclassical contributions namely, extensional 
and bending related out-of-plane warping on the accuracy of the response predictions 
is shown to be significant. Moreover, the contribution of shear deformation is shown 
to be implicitly accounted for through the bending related out-of-plane warping, and 
in-plane warping effect is found to be negligible. 

5.6 Closing Remarks 

For anisotropic beams, the major reason for the discrepancy in the predictions of the 
analytical models of Refs. [30] and [41]- [46] and the present theory is due to the apriori 
assumed displacement fields which neglect the extension and bending-related out-of- 
plane warping. The influence of the material’s anisotropy on the displacement is too 
complex to cast in a kinematic assumption similar to classical theory of extension- 
bending and torsion. 

A consistent approach to account for the various behavioral modes associated 
with anisotropic beams was adopted in this work. It is based on an asymptotical 
analysis of the energy. The influence of the material’s anisotropy on the displacement 
and stiffness coefficients was isolated, and by comparison an assessment of previous 
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analyses was performed. In particular, this approach accounts implicitly the shear 
deformation contribution shown to be significant in previous models. The difference 
being the consistent order of magnitude that this contribution is accounted for and 
its significance relative to other contributions. 
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Figure 5.1: Beam Ctoss Section 



End 


Figure 5.2: Significance of out-of-plane bending related warping on the bending slope 
of an antisymmetric [15)6 cantilever under 1 lb transverse tip Load 


Bending Slope (radians) 


no 


0.025 

0.02 


0.015 

0.01 





End 


Figure 5.3: Significance of out-of-plane bending related warping on the bending slope 
of a sj'mmetric [30J 6 cantilever under 1 lb transverse tip Load 
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Figure 5.4: Bending Slope of an Anti- Symmetric [15] 6 Cantilever Under 1 lb Trans- 
verse Tip Load 
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Figure 5.5: Bending Slope of a Symmetric 130] 6 Cantilever Under 1 lb Transverse Tip 
Load 
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Figure 5.6: Twist of a Symmetric [30]e Cantilever Under 1 lb Transverse Tip Load 
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■ Experimental [ 42 , 46] Q Smith and Chopra [42, 46] 

0 Present 13 Rehfield and Atilgan [43], Atilgan [44] 



Figure 5.9: Twist at Mid-Span Under Unit Tip Torque of Symmetric Lay-up Can- 
tilever Beams 
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Figure 5.10: Deflection of an Antisymmetric [15)6 Cantilever under 1 lb transverse 





Shear Strain x 10 4 


8 



Figure 5.11: Shear Deformation -) xy of an Antisymmetric [15] 6 Cantilever under 1 lb 
transverse tip load 
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CHAPTER VI 

CONCLUSIONS AND RECOMMENDATIONS 


This research addresses two key issues for the continuing implementation of com- 
posites in advanced structures namely, the understanding of the role of the material's 
anisotrop3 r on its stiffness behavior and its damage modes. An analytical model based 
upon a shear deformation theorj’ and a sublaminate approach was developed in or- 
der to investigate mid-plane and matrix crack-tip delaminations. This model was 
combined with an earlier analysis for mixed-mode free-edge delamination to form an 
integrated code for the prediction of damage onset in laminated composites. The 
code predictions were validated by comparing its results with test data. Of signif- 
icance is the ability it provides for the prediction of damage progression sequence 
and corresponding critical strains. Moreover, the effect of hygrothermal stresses on 
the strain energy release rate and interlaminar stresses was isolated. The increase 
in strain energy release rate and interlaminar stresses associated with curing stresses 
can precipitate failure and should be considered for an accurate prediction of failure. 

The findings of this research work point to new research inquiries. The first is 
characterization and prediction of damage onset and growth under cyclic loading 
including the effect of hygrothermal stresses. The investigation can lead to the deter- 
mination of composite components’ life and inspection intervals. The second is the 
study of the effect of damage modes and their interactions on the vibration charac- 
teristics and damping of laminated composites. The result of this investigation will 
assess the effect of damage modes on the natural frequencies and mode shapes and 
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can lead to the development of Non-Destructive Evaluation methods. 

The asymptotical analysis used to develop the thin- walled anisotropic beam theory 
provides a rigorous basis for the prediction of the beam stiffnesses and associated 
displacement field. Closed-form expressions for the stiffnesses have been developed 
and new contributions to the warping have been found. This analysis can be extended 
to beams with multi-cell type cross sections and pretwisted configurations. Moreover, 
the previous results on the effects of hygrothermal stresses point to the significance 
of including their contribution in the thin-walled closed section beam analysis. The 
consideration of dynamic and aerodynamic loadings using asymptotical analysis will 
provide a rigorous basis for the investigation of the dynamic and aeroelastic response 
of composite structures. Finally, the presence of embedded delamination on the 
response of composite beams is a first step toward studying the effect of damage 
modes on their stiffness and strength. In this respect, the analysis of composite 
beams with open cross section can be regarded as the final stage of damage in a 
closed section beam. 

When accomplished, these recommended research tasks will provide an under- 
standing of the effects of damage on the performance of advanced structures made 
out of composite and will lead to the development of reliable design tools to ensure 
their damage tolerance. 
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Appendix A 


Convergence of Displacement Field 


In this appendix detailed calculation of the third and final cycle is provided. 
Results show that no additional correction terms of the same order in the energy 
functional emerge and the displacement field given in Eq. (4.68) is the converged one. 
1.1 Third-Order Approximation 

A third cycle is carried out by rewriting the displacement field in Eq. (4.68) in the 
form 


t’j = fM*) - v(s)U! 1 (t) - r(s)f^(r) + G(s)ip'{x) 

+0l(s)f' r i'(*) + 92{s)U2(T) -f 93{s)l : 3(x) + U’i(s.x) 

t’2 = Uii*)-?- + Uaixj-j- + ( i c ’( a ') r n + w» 2 (s, a*) 
ds as 

dz dy 

v = U 3(1)3 <p(T)r t -f w(s,x) 

ds ds 


(A - 1) 


where tZ>! , tZ» 2 and w are correction functions to be determined based on their contri- 
butions to the energy functional. 

Substitute Eq. (A-l) into (4.7) to obtain the strains and curvatures in terms of 
the displacement corrections 
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where 5 0/ ? and P a & are the strains and curvatures corresponding to the 6econd-order 
approximation. These are expressed as 

ML ML ML ( ? } 

(f?> (^) (^) 

+ gAs)U"(x) + 92V? (x ) + g 3 U^'( 7 ) 

(f) mm MM MM 

25 12 = ^c^'(x)+ ^(*)+ ^r/"(x)+ 

7 22 = 0 

Pu = u ^ x )t ~ u *i x )~r - { p"( x ) Tt 0 (7^) ( A ~ 3 ) 

P22 = 0 

An order of magnitude comparison for each strain and curvature measure shows that 
some terms of higher order in 7 U can be cancelled and its expression simplifies to 

ill Jffl JjL l 1 gl 

t„ = !/;<*) - »d?(*) - zoj(*) +c(.v"(*) 

Among the newterms introduced by the function tZ>,- the leading ones are denoted 
by superscript * in Eq. (A-2). The order of tD, is assumed to be 


tli, ~ 0 


(¥) 


(A -4) 
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Consequently, the order of magnitude of the leading terms in Eq. (A-2), is as follows 


r , n( Ad2 \ 

12 ^ "1 22 ~ 0 y — j 

1 - n ( Ad \ 

P12 ~ P22 ~ 0 I -jj 1 


(A - 5 ) 


The energy functional can be represented by $ (711,2712, ->22, By 
keeping the strains and curvature associated with the second-order approximation 
and the leading terms contribution over the other terms (i.e., by dropping the terms 


§jr- and ~ TRlte in Ec b (A- 2 )) the energy function can be written as 
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w 

In the following, the order of magnitude of the energy due to bending, i.e. due to p u . 

w * ** 

p 12 , / 5 12 , and p 22 . is assessed. 

The interaction terms associated with p u , namely 
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are of order or smaller. They are neglected in comparison with the following 

membrane contribution to the energy 
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The interaction terms due to the bending curvature p J2 are 
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These terms are of higher order of magnitude in comparison with the corresponding 
membrane contribution in Eq. (A-7), and may be neglected. The remaining interac- 
tion terms associated with p 12 and p 22 , namely 


^Infill > ? ^711^22 » ^7 12^22 S 



associated with U[ and i p' 
associated with U 2 and U£ 


may also be neglected in comparison with (A-7). Therefore in order to determine the 
functions w, one has to minimize the shell energy expressed by 


I = J ^$(5 n » 2 ^i2 + 2 7i2O22>0,0,0 )dsdx (A - 8) 


Setting the first variation of the energy functional to zero to get Eq. (4.45). Sub- 
stitute from Eq. (A-2) into Eq. (4.45) to obtain 
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Equation ( A-9) shows that the contribution of w is of higher order in comparison with 
all other terms and may be cancelled from the left hand side. Therefore no additional 
corrections to the displacement field emerges, and the displacement field obtained in 
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Eq. (4.68) is the converged one. An alternative is to neglect the terms of higher order 
in Eq. (A-9). while keeping the leading u>i term, to obtain 






2A e , dg\ T , t N dg? „ dgz „ dii'i 

— (x) 4 — (r,(x) 4 — t/ s (x) + _£7 s(l ) + _ 


= constant 


(A -10) 

Solution of Eq. (A- 10) is determined using the single value condition of the axial 
displacement and Wi is found to be a function of x onl}-. Such a function has already 
been considered and no new terms of the same order in the energy functional are 
generated from the third and therefore final cycle. 
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